# %%writefile requirements.txt
# pandas
# numpy
# matplotlib
# geopandas
# cartopy
# rasterio
# atlite
# networkx
# pypsa
# highspy
# xarray
# rioxarray

# from google.colab import drive
# drive.mount('/content/drive')

# !pip install -r requirements.txt

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
#import cartopy 
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import rasterio as rio 
import atlite
from atlite.gis import ExclusionContainer
from atlite.gis import shape_availability
from rasterio.plot import show
import matplotlib.colors as mcolors
import matplotlib.patches as mpatches
import networkx as nx
import pypsa
import country_converter as coco
import xarray as xr
import os

path = "gadm_410-levels-ADM_1-NLD.gpkg"

country = gpd.read_file(path)

country.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

fig, ax = plt.subplots(subplot_kw={'projection': ccrs.EqualEarth()}, figsize=(10, 10))

country.plot(ax=ax, edgecolor="black", facecolor = "lightgrey", linewidth = 0.2)

<GeoAxes: >

ocean = gpd.read_file("eez_v11.gpkg")

ocean.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

ocean_nld = ocean[ocean["SOVEREIGN1"] == "Netherlands"]
ocean_nld

ocean_nld["GEONAME"]

61             Aruban Exclusive Economic Zone
66            Bonaire Exclusive Economic Zone
74          Curaçaoan Exclusive Economic Zone
92               Saba Exclusive Economic Zone
94     Sint-Eustatius Exclusive Economic Zone
100      Sint-Maarten Exclusive Economic Zone
216             Dutch Exclusive Economic Zone
Name: GEONAME, dtype: object

fig, ax = plt.subplots(subplot_kw={'projection': ccrs.EqualEarth()}, figsize=(10, 10))

ocean_nld[ocean_nld.index.isin([216])].plot(ax=ax, edgecolor="black", linewidth = 0.2)

<GeoAxes: >

eez_nld_4326 = ocean_nld.copy()
eez_nld_4326 = eez_nld_4326[eez_nld_4326["GEONAME"] == "Dutch Exclusive Economic Zone"]
eez_nld_4326["region"] = [6]

fig, ax = plt.subplots(subplot_kw={'projection': ccrs.EqualEarth()})

country.plot(ax=ax, edgecolor="black", facecolor="lightgrey",linewidth = 0.2)
ocean_nld[ocean_nld.index.isin([216])].plot(ax=ax, edgecolor="black", facecolor="lightblue",linewidth = 0.2)

<GeoAxes: >

nld_4326 = country.copy()
nld_4326

fig, ax = plt.subplots()
nld_4326.plot(ax=ax)

#Visualing the different regions of spain with their index
for idx, geometry in enumerate(nld_4326.geometry):
    centroid = geometry.centroid
    ax.annotate(str(idx), (centroid.x, centroid.y))

plt.show()

nld_4326["region"] = [1,2,1,2,1,6,3,3,4,2,4,5,6,5]
nld_4326_dis = nld_4326.dissolve(by="region")
nld_4326_dis

rep_point = nld_4326_dis.representative_point()
colors =  ('orange', 'skyblue', 'lightgreen', 'purple', 'cyan', 'blue')

fig, ax = plt.subplots()
nld_4326_dis.plot(ax=ax, color = colors, edgecolor="black", linewidth=0.2)  # Tracer les géométries

# Ajouter les numéros d'index comme annotations
for idx, geometry in enumerate(nld_4326_dis.geometry):
    centroid = geometry.centroid
    ax.annotate(str(idx+1), (centroid.x, centroid.y))

plt.show()

merged_df = pd.merge(nld_4326_dis, eez_nld_4326, on=['region', 'geometry'], how='outer')
merged_df_dis = merged_df.dissolve(by="region")

merged_df_dis

# Définir les couleurs pour chaque région
regions = merged_df_dis.index.unique()
num_regions = len(regions)
cmap = plt.get_cmap('tab20')
colors = [cmap(i) for i in np.linspace(0, 1, num_regions)]
color_dict = dict(zip(regions, colors))
rep_point = merged_df_dis.representative_point()

# Plot
fig, ax = plt.subplots(1, 1, figsize=(10, 6))

for region in regions:
    subset = merged_df_dis[merged_df_dis.index == region]
    subset.plot(ax=ax, color=color_dict[region], label=region, edgecolor = "black", linewidth = 0.2)

rep_point.plot(ax=ax, color='red', markersize=3)
for idx, region in enumerate(merged_df_dis.index):
    ax.annotate(region, (rep_point.x[region], rep_point.y[region]))

# Ajouter la légende
handles = [plt.Line2D([0], [0], marker='s', color=color_dict[region], linestyle='', markersize=10, label=region) for region in regions]
ax.legend(handles=handles, title='Region', bbox_to_anchor=(1.05, 1), loc='upper left')

plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('GeoDataFrame plot by region')
plt.show()

nld_land_4326 = merged_df_dis.loc[[1,2,3,4,5]]
nld_sea_4326 = merged_df_dis.loc[[6]]

url_elevation = "GEBCO_2014_2D-NL.nc"
url_airports = "ne_10m_airports.gpkg"
url_roads = "ne_10m_roads.gpkg"
url_nature = "WDPA_Oct2022_Public_shp-NLD.tif"
url_land = "PROBAV_LC100_global_v3.0.1_2019-nrt_Discrete-Classification-map_EPSG-4326-NL.tif"
fn = "gadm_410-levels-ADM_1-NLD.gpkg"

countries_on = nld_land_4326
shape_on = countries_on.to_crs(epsg=3035).geometry

shape_on.plot()

<Axes: >

excluder_onwind = ExclusionContainer(crs=3035, res = 300)

# Protected areas excluded
excluder_onwind.add_raster(url_nature, crs=3035)

# Airports excluded
excluder_onwind.add_geometry(url_airports, buffer=10000)

# Roads excluded
excluder_onwind.add_geometry(url_roads, buffer = 300)

# Elevation above 2000 excluded
excluder_onwind.add_raster(url_elevation, codes=lambda x: x > 2000, crs=4326)

# Urban areas excluded (see Copernicus classification above)
excluder_onwind.add_raster(url_land, codes=[50], buffer=1000, crs=4326)

# Areas to include
excluder_onwind.add_raster(url_land, codes=[20,30,40,60,125], crs=4326, invert=True)

band, transform = shape_availability(shape_on, excluder_onwind)

# geometry_0 = shape_on.iloc[0]

# # Créer un masque pour la région 0
# mask = rio.features.geometry_mask([geometry_0], out_shape=band.shape, transform=transform, invert=True)

# # Compter le nombre de cases True dans band correspondant à la région 0
# num_true_in_region_0 = np.sum(band[mask])

# shape_on

# for index, (key, value) in enumerate(shape_on.items()):
#     print(f"Index: {index}, Key: {key}, Value: {value}")

density_installation = {
    'onwind': 5,    # potential MW/km2
    'offwind': 7, 
    'solar': 100    
}

def renewable_potential(mask, tr, excl, shp, tech) : 

    if tech not in density_installation:
        raise ValueError(f"Unknown technology : {tech}")
    else : 

        potentials = []
        keys = []
        density = density_installation[tech]
        cell_area = excl.res**2

        for i, (region, value) in enumerate(shp.items()):  
            geom = shp.iloc[i]
            mask_r = rio.features.geometry_mask([geom], out_shape=mask.shape, transform=tr, invert=True)
            el_cel_in_r = np.sum(mask[mask_r])
            el_area_in_r = el_cel_in_r * cell_area
            el_area_in_r_km2 = el_area_in_r / 1e6
            potential = el_area_in_r_km2 * density
            potentials.append(potential)
            keys.append(region)
        
        dict_pot = {keys[i]: potentials[i] for i in range(len(keys))}
        
        return dict_pot

renewable_potential(band, transform, excluder_onwind, shape_on, "onwind")

{1: 11943.45,
 2: 10357.650000000001,
 3: 3581.1000000000004,
 4: 2205.45,
 5: 3563.1}

fig, ax = plt.subplots(figsize=(10, 6))
shape_on.plot(ax=ax, color="none", linewidth = 0.2)
show(band, transform=transform, cmap="Greens", ax=ax)

eligible_patch = mpatches.Patch(facecolor='green', label='Eligible Areas')
non_eligible_patch = mpatches.Patch(facecolor='white', label='Non-Eligible Areas', edgecolor='black')

plt.legend(handles=[eligible_patch, non_eligible_patch], loc='upper right')

plt.title("Eligible Areas For Wind Turbins")
plt.show()

# ocean = gpd.read_file("/content/drive/MyDrive/DataDSEM/eez_v11.gpkg")
# ocean_nl = ocean[ocean["SOVEREIGN1"] == "Netherlands"]
# ocean_nl.crs
# ocean_bound = gpd.read_file("/content/drive/MyDrive/DataDSEM/eez_boundaries_v11.gpkg")
# bound_esp = ocean_bound[ocean_bound["SOVEREIGN1"] == "Spain"]

countries_off = nld_sea_4326
shape_off = countries_off.to_crs(epsg=3035).geometry

shape_off.plot()

<Axes: >

excluder_offshore = ExclusionContainer(crs=3035, res = 300)

# Exclude protected areas
excluder_offshore.add_raster(url_nature, crs=3035)

# Exclude depth beneath 50m
excluder_offshore.add_raster(url_elevation, codes=lambda x: x<(-50), crs=4326)

# Away from coutry shores
excluder_offshore.add_geometry(nld_4326_dis, buffer=10000)

excluder_offshore.add_geometry('eez_boundaries_v11.gpkg') # Within EEZ boundaries

band_off, transform_off = shape_availability(shape_off, excluder_offshore)

fig, ax = plt.subplots(figsize=(10, 6))
shape_off.plot(ax=ax, color="none", linewidth = 0.2)
show(band_off, transform=transform_off, cmap="Greens", ax=ax)

eligible_patch = mpatches.Patch(facecolor='green', label='Eligible Areas')
non_eligible_patch = mpatches.Patch(facecolor='white', label='Non-Eligible Areas', edgecolor='black')

plt.legend(handles=[eligible_patch, non_eligible_patch], loc='upper right')

plt.title("Eligible Areas For Wind Turbins")
plt.show()

countries_so = nld_land_4326
shape_so = countries_so.to_crs(epsg=3035).geometry

excluder_solar = ExclusionContainer(crs=3035, res = 300)

# Exclude protected areas
excluder_solar.add_raster(url_nature, crs=3035)

# Exclude land areas
excluder_solar.add_raster(url_land, codes=[50,70,80,90]+[i for i in range(111,127)], buffer=300, crs=4326)

# Eigible land areas to include
excluder_solar.add_raster(url_land, codes=[20,30,40,60,100], crs=4326, invert=True)

band_so, transform_so = shape_availability(shape_so, excluder_solar)

fig, ax = plt.subplots(figsize=(10, 6))
shape_so.plot(ax=ax, color="none", linewidth = 0.2)
show(band, transform=transform, cmap="Greens", ax=ax)

eligible_patch = mpatches.Patch(facecolor='green', label='Eligible Areas')
non_eligible_patch = mpatches.Patch(facecolor='white', label='Non-Eligible Areas', edgecolor='black')

plt.legend(handles=[eligible_patch, non_eligible_patch], loc='upper right')

plt.title("Non Suitable Areas For Solar Farms")
plt.show()

# Calculate the maximum potentials 

pot_max_onwind = renewable_potential(band, transform, excluder_onwind, shape_on, "onwind")
pot_max_offwind = renewable_potential(band_off, transform_off, excluder_offshore, shape_off, "offwind")
pot_max_solar = renewable_potential(band_so, transform_so, excluder_solar, shape_so, "solar")

print(pot_max_onwind, pot_max_offwind, pot_max_solar)

{1: 11943.45, 2: 10357.650000000001, 3: 3581.1000000000004, 4: 2205.45, 5: 3563.1} {6: 316832.04000000004} {1: 233298.0, 2: 131814.0, 3: 42201.0, 4: 61389.0, 5: 101124.0}

from urllib.request import urlretrieve

# fn = "era5-2013-NL.nc"
# url = "https://tubcloud.tu-berlin.de/s/bAJj9xmN5ZLZQZJ/download/" + fn
# urlretrieve(url, fn)

# fn = "gadm_410-levels-ADM_1-NLD.gpkg"
# regions = gpd.read_file(fn)
# regions.head()

# regions.crs

# regions.total_bounds

# nld = regions.copy()
# nld["region"] = [1,2,1,2,1,2,3,3,2,1,2,3,3,3]
# nld_dis = nld.dissolve(by = "region")
# minx, miny, maxx, maxy = nld_dis.total_bounds
# buffer = 0.25

# nld_dis.total_bounds

minx, miny, maxx, maxy = nld_land_4326.total_bounds
buffer = 0.25

cutout = atlite.Cutout("europe-2013-era5.nc")
cutout_nld = cutout.sel(x=slice(minx - buffer, maxx + buffer), y=slice(miny - buffer, maxy + buffer))

cutout_nld

<Cutout "europe-2013-era5-72w3w1c0">
 x = 3.30 ⟷ 7.20, dx = 0.30
 y = 50.70 ⟷ 53.70, dy = 0.30
 time = 2013-01-01 ⟷ 2013-12-31, dt = h
 module = era5
 prepared_features = ['height', 'wind', 'influx', 'temperature', 'runoff']

cutout_nld.data

<xarray.Dataset> Size: 86MB
Dimensions:           (x: 14, y: 11, time: 8760)
Coordinates:
  * x                 (x) float64 112B 3.3 3.6 3.9 4.2 4.5 ... 6.3 6.6 6.9 7.2
  * y                 (y) float64 88B 50.7 51.0 51.3 51.6 ... 53.1 53.4 53.7
  * time              (time) datetime64[ns] 70kB 2013-01-01 ... 2013-12-31T23...
    lon               (x) float64 112B dask.array<chunksize=(14,), meta=np.ndarray>
    lat               (y) float64 88B dask.array<chunksize=(11,), meta=np.ndarray>
Data variables: (12/13)
    height            (y, x) float32 616B dask.array<chunksize=(11, 14), meta=np.ndarray>
    wnd100m           (time, y, x) float32 5MB dask.array<chunksize=(100, 11, 14), meta=np.ndarray>
    wnd_azimuth       (time, y, x) float32 5MB dask.array<chunksize=(100, 11, 14), meta=np.ndarray>
    roughness         (time, y, x) float32 5MB dask.array<chunksize=(100, 11, 14), meta=np.ndarray>
    influx_toa        (time, y, x) float32 5MB dask.array<chunksize=(100, 11, 14), meta=np.ndarray>
    influx_direct     (time, y, x) float32 5MB dask.array<chunksize=(100, 11, 14), meta=np.ndarray>
    ...                ...
    albedo            (time, y, x) float32 5MB dask.array<chunksize=(100, 11, 14), meta=np.ndarray>
    solar_altitude    (time, y, x) float64 11MB dask.array<chunksize=(100, 11, 14), meta=np.ndarray>
    solar_azimuth     (time, y, x) float64 11MB dask.array<chunksize=(100, 11, 14), meta=np.ndarray>
    temperature       (time, y, x) float64 11MB dask.array<chunksize=(100, 11, 14), meta=np.ndarray>
    soil temperature  (time, y, x) float64 11MB dask.array<chunksize=(100, 11, 14), meta=np.ndarray>
    runoff            (time, y, x) float32 5MB dask.array<chunksize=(100, 11, 14), meta=np.ndarray>
Attributes:
    module:             era5
    prepared_features:  ['height', 'runoff', 'temperature', 'wind', 'influx']
    chunksize_time:     100
    dx:                 0.3
    dy:                 0.3
    Conventions:        CF-1.6
    history:            2022-03-04 14:08:37 GMT by grib_to_netcdf-2.24.2: /op...

array([3.3, 3.6, 3.9, 4.2, 4.5, 4.8, 5.1, 5.4, 5.7, 6. , 6.3, 6.6, 6.9, 7.2])

array([50.7, 51. , 51.3, 51.6, 51.9, 52.2, 52.5, 52.8, 53.1, 53.4, 53.7])

array(['2013-01-01T00:00:00.000000000', '2013-01-01T01:00:00.000000000',
       '2013-01-01T02:00:00.000000000', ..., '2013-12-31T21:00:00.000000000',
       '2013-12-31T22:00:00.000000000', '2013-12-31T23:00:00.000000000'],
      dtype='datetime64[ns]')

PandasIndex(Index([3.3, 3.6, 3.9, 4.2, 4.5, 4.8, 5.1, 5.4, 5.7, 6.0, 6.3, 6.6, 6.9, 7.2], dtype='float64', name='x'))

PandasIndex(Index([50.7, 51.0, 51.3, 51.6, 51.9, 52.2, 52.5, 52.8, 53.1, 53.4, 53.7], dtype='float64', name='y'))

PandasIndex(DatetimeIndex(['2013-01-01 00:00:00', '2013-01-01 01:00:00',
               '2013-01-01 02:00:00', '2013-01-01 03:00:00',
               '2013-01-01 04:00:00', '2013-01-01 05:00:00',
               '2013-01-01 06:00:00', '2013-01-01 07:00:00',
               '2013-01-01 08:00:00', '2013-01-01 09:00:00',
               ...
               '2013-12-31 14:00:00', '2013-12-31 15:00:00',
               '2013-12-31 16:00:00', '2013-12-31 17:00:00',
               '2013-12-31 18:00:00', '2013-12-31 19:00:00',
               '2013-12-31 20:00:00', '2013-12-31 21:00:00',
               '2013-12-31 22:00:00', '2013-12-31 23:00:00'],
              dtype='datetime64[ns]', name='time', length=8760, freq=None))

#nld_dis.index = nld_dis.GID_1

#nld_land_4326.crs

onshore_A = cutout_nld.availabilitymatrix(nld_land_4326, excluder_onwind)

Compute availability matrix:   0%|          | 0/5 [00:00<?, ? gridcells/s]

Compute availability matrix: 100%|██████████| 5/5 [01:16<00:00, 15.37s/ gridcells]

cap_per_sqkm = 3  # MW/km2

area = cutout_nld.grid.set_index(["y", "x"]).to_crs(3035).area / 1e6

area = xr.DataArray(area, dims=("spatial"))

capacity_matrix = onshore_A.stack(spatial=["y", "x"]) * area * cap_per_sqkm

onshore = cutout_nld.wind(
    atlite.windturbines.Vestas_V112_3MW,
    matrix=capacity_matrix,
    index=nld_land_4326.index,
    per_unit=True,
)

onshore_df = onshore.to_pandas()
onshore_df

/opt/anaconda3/envs/pitie/lib/python3.11/site-packages/atlite/resource.py:72: FutureWarning: 'add_cutout_windspeed' for wind turbine
power curves will default to True in atlite relase v0.2.13.
  warnings.warn(msg, FutureWarning)
INFO:atlite.convert:Convert and aggregate 'wind'.

[########################################] | 100% Completed | 4.38 s

region_names = nld_land_4326.index
for region in region_names:
   
    us = pd.Series(onshore_df[region], index=onshore_df.index)
    

    fig, ax = plt.subplots(figsize=(15,5))
    plt.title(f'Capacity factors in {region} for the month of January 2013')
    us.loc['2013-01'].plot(color="blue", linestyle='-.', label='Onshore')
    plt.ylabel('Capacity Factor')
    plt.legend(loc="upper right")

minx, miny, maxx, maxy = nld_sea_4326.total_bounds
buffer = 0.25

cutout = atlite.Cutout("europe-2013-era5.nc")
cutout_ocean = cutout.sel(x=slice(minx - buffer, maxx + buffer), y=slice(miny - buffer, maxy + buffer))

offshore_A = cutout_ocean.availabilitymatrix(nld_sea_4326, excluder_offshore)

Compute availability matrix: 100%|██████████| 1/1 [00:01<00:00,  1.67s/ gridcells]

cap_per_sqkm = 3  # MW/km2

area = cutout_ocean.grid.set_index(["y", "x"]).to_crs(3035).area / 1e6

area = xr.DataArray(area, dims=("spatial"))

capacity_matrix = offshore_A.stack(spatial=["y", "x"]) * area * cap_per_sqkm

#create offshore wind profiles
offshore = cutout_ocean.wind(
    atlite.windturbines.NREL_ReferenceTurbine_5MW_offshore,
    matrix= capacity_matrix,
    index= nld_sea_4326.index,
    per_unit=True,
)

offshore_df = offshore.to_pandas()

/opt/anaconda3/envs/pitie/lib/python3.11/site-packages/atlite/resource.py:72: FutureWarning: 'add_cutout_windspeed' for wind turbine
power curves will default to True in atlite relase v0.2.13.
  warnings.warn(msg, FutureWarning)
INFO:atlite.convert:Convert and aggregate 'wind'.

[########################################] | 100% Completed | 27.57 s

offshore_df

region_names = nld_sea_4326.index

for region in region_names:
   
    vs = pd.Series(offshore_df[region], index=offshore_df.index)
    
    fig, ax = plt.subplots(figsize=(15,5))
    plt.title(f'Capacity factors in {region} for the month of January 2013')
    us.loc['2013-01'].plot(color="green", linestyle='-.', label='Offshore')
    plt.ylabel('Capacity Factor')
    plt.legend(loc="upper right")

minx, miny, maxx, maxy = nld_land_4326.total_bounds
buffer = 0.25

cutout = atlite.Cutout("europe-2013-era5.nc")
cutout_nld = cutout.sel(x=slice(minx - buffer, maxx + buffer), y=slice(miny - buffer, maxy + buffer))

solar_A = cutout_nld.availabilitymatrix(nld_land_4326, excluder_solar)

Compute availability matrix: 100%|██████████| 5/5 [00:01<00:00,  2.56 gridcells/s]

cap_per_sqkm = 3  # MW/km2

area = cutout_nld.grid.set_index(["y", "x"]).to_crs(3035).area / 1e6

area = xr.DataArray(area, dims=("spatial"))

capacity_matrix = solar_A.stack(spatial=["y", "x"]) * area * cap_per_sqkm

#create solar profile
solar = cutout_nld.pv(
    panel=atlite.solarpanels.CdTe,
    matrix=capacity_matrix,
    orientation='latitude_optimal',
    index= nld_land_4326.index,
    per_unit=True,
)

solar_df = solar.to_pandas()
solar_df

INFO:atlite.convert:Convert and aggregate 'pv'.

[########################################] | 100% Completed | 12.08 s

solar_df[2]

time
2013-01-01 00:00:00    0.0
2013-01-01 01:00:00    0.0
2013-01-01 02:00:00    0.0
2013-01-01 03:00:00    0.0
2013-01-01 04:00:00    0.0
                      ... 
2013-12-31 19:00:00    0.0
2013-12-31 20:00:00    0.0
2013-12-31 21:00:00    0.0
2013-12-31 22:00:00    0.0
2013-12-31 23:00:00    0.0
Name: 2, Length: 8760, dtype: float64

region_names = nld_land_4326.index

for region in region_names:
   
    us = pd.Series(solar_df[region], index=solar_df.index)
    

    fig, ax = plt.subplots(figsize=(15,5))
    plt.title(f'Capacity factors in {region} for the month of January 2013')
    us.loc['2013-01'].plot(color="red", linestyle='-.', label='Solar')
    plt.ylabel('Capacity Factor')
    plt.legend(loc="upper right")

df_offwind_cf = offshore_df.copy()
df_onwind_cf = onshore_df.copy()
df_solar_cf = solar_df.copy()

df_offwind_cf = df_offwind_cf.drop(labels = df_offwind_cf.loc['2020-02-29 00:00:00':'2020-02-29 23:00:00'].index.array , axis = 0)
df_onwind_cf = df_onwind_cf.drop(labels = df_onwind_cf.loc['2020-02-29 00:00:00':'2020-02-29 23:00:00'].index.array , axis = 0)
df_solar_cf = df_solar_cf.drop(labels = df_solar_cf.loc['2020-02-29 00:00:00':'2020-02-29 23:00:00'].index.array , axis = 0)
len(df_solar_cf) == len(df_onwind_cf) == len(df_offwind_cf) == 8760

True

df_pp = pd.read_csv("global_power_plant_database.csv")
gdf_pp = gpd.GeoDataFrame(df_pp, geometry=gpd.points_from_xy(df_pp.longitude, df_pp.latitude))
gdf_pp.set_crs(epsg=4326, inplace=True)

#only NLP pp
gdf_pp_nl = gdf_pp[gdf_pp.country == 'NLD']

#No existing wind and solar pp's to be considered
#gdf_pp_nl = gdf_pp_nl.drop(gdf_pp_nl[(gdf_pp_nl.primary_fuel == "Solar")].index)
#gdf_pp_nl = gdf_pp_nl.drop(gdf_pp_nl[(gdf_pp_nl.primary_fuel == "Wind")].index)

/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/ipykernel_27816/2596194382.py:1: DtypeWarning: Columns (10) have mixed types. Specify dtype option on import or set low_memory=False.
  df_pp = pd.read_csv("global_power_plant_database.csv")

gdf_pp_nl.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

gdf_pp_nl.head()

unique_values = gdf_pp_nl['primary_fuel'].unique()
unique_values

array(['Solar', 'Wind', 'Nuclear', 'Gas', 'Coal'], dtype=object)

dff = merged_df_dis.copy()
dff = dff.reset_index()
dff

# df_nld = nld_land_4326.reset_index()
# df_onld = nld_sea_4326.reset_index()
# merged_df = pd.merge(df_nld, df_onld, on=['region', 'geometry'], how='outer')

# Join in such a way that the geometry of the powerplant is included in the geomtetry of the region 
joined = gpd.sjoin(gdf_pp_nl, dff, how='left', predicate='within')
joined

joined.head()
joined.to_csv('pp_nld_region.csv')

df_fossils = joined.groupby(["region","primary_fuel"]).capacity_mw.sum().unstack().fillna(0)
df_fossils

cap_pp_r = df_fossils.copy()
cap_pp_r.index = df_fossils.index.astype(int)
cap_pp_r

merged = merged_df_dis.join(cap_pp_r, how='inner')
merged

merged.representative_point()

region
1    POINT (6.30770 53.03945)
2    POINT (6.29240 52.29567)
3    POINT (5.83110 51.27609)
4    POINT (4.80150 52.41195)
5    POINT (4.41346 51.97036)
6    POINT (4.13744 53.40210)
dtype: geometry

# plt.figure(figsize=(8, 8))
# plt.pie(merged['Coal'], labels=merged.index, autopct='%1.1f%%', startangle=140, colors=['blue', 'orange', 'green', 'red', 'purple', 'pink'])
# plt.title('Capacities by Type of Energy')
# plt.show()

merged.representative_point().iloc[1].coords[0]

(6.292400568002548, 52.29567146350013)

merged['geometry'].iloc[1].representative_point().coords[0]

(6.292400568002548, 52.29567146350013)

bar_plot = merged[['Coal','Gas','Nuclear']].hvplot.bar(
    stacked=True,  # Stack the bars to show distribution within each region
    xlabel='Region',
    ylabel='Capacity (MW)',
    title='Distribution of Energy Capacities by Region',
    rot=45,  # Rotate x-axis labels for better readability
    legend='top_right'
)

bar_plot.opts(legend_position='top_right', title="Distribution of Existing Powerplants' Capacities by Region")
bar_plot

# data = {
#     'Region': [1, 2, 3],
#     'Coal': [100, 150, 200],
#     'Gas': [200, 250, 300],
#     'Nuclear': [300, 350, 400]
# }
# df = pd.DataFrame(data)

# # Utilisation de iterrows et row.get
# for idx, row in df.iterrows():
#     #print(f"Index: {idx}")
#     #print(f"Row data:\n{row}")
    
#     coal_capacity = row.get('Coal')
#     gas_capacity = row.get('Gas')
#     nuclear_capacity = row.get('Nuclear')
#     non_existent = row.get('NonExistentColumn', 'DefaultValue')
    
#     print(f"Coal capacity: {coal_capacity}")
#     print(f"Gas capacity: {gas_capacity}")
#     print(f"Nuclear capacity: {nuclear_capacity}")
#     print(f"Non-existent column value: {non_existent}")
#     print("-" * 40)

# merged.to_csv('testte.csv')

# gdf = merged.copy()
# gdf['representative_point'] = gdf.geometry.representative_point()

# # Créer la carte
# fig, ax = plt.subplots(figsize=(12, 12))
# gdf.plot(ax=ax, color='blue', edgecolor='black', label='Geometries')
# gdf.set_geometry('representative_point').plot(ax=ax, color='red', markersize=5, label='Representative Points')

# # Ajouter une légende
# plt.legend()

# # Afficher la carte
# plt.title('Geometries and Their Representative Points')
# plt.xlabel('Longitude')
# plt.ylabel('Latitude')
# plt.show()

import holoviews as hv

import hvplot.pandas

df_pp = pd.read_csv("global_power_plant_database.csv")
gdf_pp = gpd.GeoDataFrame(df_pp, geometry=gpd.points_from_xy(df_pp.longitude, df_pp.latitude))
gdf_pp.set_crs(epsg=4326, inplace=True)

/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/ipykernel_27816/3143510209.py:1: DtypeWarning: Columns (10) have mixed types. Specify dtype option on import or set low_memory=False.
  df_pp = pd.read_csv("global_power_plant_database.csv")

gdf_pp_nl = gdf_pp[gdf_pp.country == 'NLD']
gdf_pp_nl

gdf_pp_nl.query("primary_fuel == 'Gas'").hvplot(
    geo=True, tiles=True
)

plot = gdf_pp_nl.query("primary_fuel in ['Nuclear', 'Coal', 'Gas']").hvplot(
    geo=True,
    tiles="CartoLight",
    frame_height=600,
    color="primary_fuel",
    cmap="Category10",
    size="capacity_mw",
    alpha=0.6,
    title="Gas, Coal and Nucear Powerplants in NLD"
)
plot

country

country[~country.index.isin([5])].plot()

<Axes: >

country["NAME_1"]

0           Drenthe
1         Flevoland
2           Fryslân
3        Gelderland
4         Groningen
5        IJsselmeer
6           Limburg
7     Noord-Brabant
8     Noord-Holland
9        Overijssel
10          Utrecht
11          Zeeland
12    Zeeuwse meren
13               NA
Name: NAME_1, dtype: object

country.plot()

<Axes: >

""""
Websites for data : 

https://www.statista.com/statistics/753196/total-number-of-inhabitants-in-the-netherlands-by-province/

https://www.cbs.nl
"""

population_nld = {
    "Drenthe" : 502051,
    "Flevoland" :444701,
    "Fryslân" : 659551,
    "Gelderland": 2133708,
    "Groningen":596075,
    "IJsselmeer":0,
    "Limburg":1128367,
    "Noord-Brabant":2626210,
    "Noord-Holland":2952622,
    "Overijssel":1184333,
    "Utrecht":1387643,
    "Zeeland": 391124,
    "Zeeuwse meren": 0,
    "NA": 3804906
    }

# population_eez = {
#     "Aruban Exclusive Economic Zone": 107609,
#     "Bonaire Exclusive Economic Zone": 22573,
#     "Curaçaoan Exclusive Economic Zone": 152369,
#     "Saba Exclusive Economic Zone": 1933,
#     "Sint-Eustatius Exclusive Economic Zone": 3242,
#     "Sint-Maarten Exclusive Economic Zone": 41486,
#     "Dutch Exclusive Economic Zone": 0
#     }

population_eez = {
    "Dutch Exclusive Economic Zone": 0
    }

df_population_eez = pd.DataFrame.from_dict(population_eez,orient="index")
df_population_eez.rename(columns = {0:'population'}, inplace = True)
df_population_eez["region"] = [6]

    
df_population_nld = pd.DataFrame.from_dict(population_nld,orient="index")
df_population_nld.rename(columns = {0:'population'}, inplace = True)
df_population_nld["region"] = [1,2,1,2,1,6,3,3,4,2,4,5,6,5]


df_population = pd.concat([df_population_nld, df_population_eez])
df_population["ratio"] = df_population.population/df_population.population.sum()
df_population

region_population = df_population.groupby("region").ratio.sum().to_frame()
region_population["population"] = df_population.groupby("region").population.sum()
region_population

df_reg = region_population['population']
df_reg

region
1    1757677
2    3762742
3    3754577
4    4340265
5    4196030
6          0
Name: population, dtype: int64

def func(pct, allvals):
    return "{:.2f}%".format(pct)

# Create the pie chart
fig, ax = plt.subplots(figsize=(6, 6))
wedges, texts, autotexts = ax.pie(df_reg, labels=df_reg.index, autopct=lambda pct: func(pct, df_reg), startangle=140, textprops=dict(color="w"))

# Format the text on the pie chart
for text in autotexts:
    text.set_color('black')
    text.set_weight('bold')

# Create custom legend with absolute values
legend_labels = [f"{label}: {value:.0f} inhabitants" for label, value in zip(df_reg.index, df_reg)]
ax.legend(wedges, legend_labels, title="Region", loc="center left", bbox_to_anchor=(1, 0, 0.5, 1))

ax.set_title("Population by region")


plt.show()

load = pd.read_csv("load.csv",index_col="time")
load = load["NL"]
load = load.to_frame()
load.head(3)
load.index = pd.to_datetime(load.index)
load

fig, axs = plt.subplots(figsize=(15, 10))

axs.bar(load.index, load['NL']/1e3)
axs.set_title('electricity demand Netherlands 2013 [GW]')
axs.set_xlabel('Day')
axs.set_ylabel('demand [GW]')

Text(0, 0.5, 'demand [GW]')

for idx in region_population.index : 
    load[idx] = load["NL"].values * region_population.loc[idx, 'ratio']

load.head(3)

regions = [i for i in range(1,7)]
mean_values = load[regions].mean()

plt.figure(figsize=(10, 6))
mean_values.plot(kind='bar', color=['blue', 'orange', 'green', 'red', 'purple', 'brown'])
plt.xlabel('Region')
plt.ylabel('Average Load')
plt.title('Average Hourly Load by Region [MW]')
plt.legend(['Region 1', 'Region 2', 'Region 3', 'Region 4', 'Region 5', 'Region 6'], title='Regions')
plt.xticks(rotation=0)
plt.tight_layout()

load.index = pd.to_datetime(load.index)
load.index

DatetimeIndex(['2013-01-01 00:00:00', '2013-01-01 01:00:00',
               '2013-01-01 02:00:00', '2013-01-01 03:00:00',
               '2013-01-01 04:00:00', '2013-01-01 05:00:00',
               '2013-01-01 06:00:00', '2013-01-01 07:00:00',
               '2013-01-01 08:00:00', '2013-01-01 09:00:00',
               ...
               '2013-12-31 14:00:00', '2013-12-31 15:00:00',
               '2013-12-31 16:00:00', '2013-12-31 17:00:00',
               '2013-12-31 18:00:00', '2013-12-31 19:00:00',
               '2013-12-31 20:00:00', '2013-12-31 21:00:00',
               '2013-12-31 22:00:00', '2013-12-31 23:00:00'],
              dtype='datetime64[ns]', name='time', length=8760, freq=None)

res = 3
load_res = load.resample(f"{res}h").first()
df_onwind_cf_res = df_onwind_cf.resample(f"{res}h").first()
df_offwind_cf_res = df_offwind_cf.resample(f"{res}h").first()
df_solar_cf_res = df_solar_cf.resample(f"{res}h").first()
df_offwind_cf_res

load_h = pd.read_csv("load.csv",index_col="time")
load_h = load_h["NL"]
load_h = load_h.to_frame()
load_h.head(3)
load_h.index = pd.to_datetime(load_h.index)

for idx in region_population.index : 
    load_h[idx] = load_h["NL"].values * region_population.loc[idx, 'ratio']

selected_columns = load_h[[1, 2, 3, 4, 5, 6]]
daily_mean = selected_columns.resample('D').mean()

# Plot the selected columns
fig, ax = plt.subplots(figsize=(12, 6))
daily_mean.plot(ax=ax)
ax.set_title('Daily load[MW] per region ')
ax.set_xlabel('Date')
ax.set_ylabel('Values')
plt.legend(title='Region', loc='lower right')
plt.show()

# Importing technology costs for 2025
year = 2025
fn = f"https://raw.githubusercontent.com/PyPSA/technology-data/master/outputs/costs_{year}.csv"
costs_df = pd.read_csv(fn, index_col=[0,1])

costs_df

# Annuity function
def annuity(r, n):
    return r / (1.0 - 1.0 / (1.0 + r) ** n)

costs = costs_df.copy()
costs.loc[costs.unit.str.contains("/kW"), "value"] *= 1e3
costs.unit = costs.unit.str.replace("/kW", "/MW")

defaults = {
    "FOM": 0,
    "VOM": 0,
    "efficiency": 1,
    "fuel": 0,
    "investment": 0,
    "lifetime": 25,
    "CO2 intensity": 0,
    "discount rate": 0.07,
}

costs = costs.value.unstack().fillna(defaults)

#Fixing costs for OCGT/CCGT
costs.at["OCGT", "fuel"] = costs.at["gas", "fuel"]
costs.at["CCGT", "fuel"] = costs.at["gas", "fuel"]
costs.at["OCGT", "CO2 intensity"] = costs.at["gas", "CO2 intensity"]
costs.at["CCGT", "CO2 intensity"] = costs.at["gas", "CO2 intensity"]

costs

annuity_value = costs.apply(lambda x: annuity(x["discount rate"], x["lifetime"]), axis=1)
costs["marginal_cost"] = costs["VOM"] + costs["fuel"] / costs["efficiency"]
costs["capital_cost"] = (annuity_value + costs["FOM"] / 100) * costs["investment"]

n = pypsa.Network()

n.snapshot_weightings.loc[:,:] = res

n.set_snapshots(load.index)

rep_points = merged_df_dis.representative_point()
rep_points_gdf = gpd.GeoDataFrame(geometry=rep_points)
rep_points_gdf['x'] = rep_points_gdf['geometry'].apply(lambda point: point.x)
rep_points_gdf['y'] = rep_points_gdf['geometry'].apply(lambda point: point.y)
rep_points_gdf

for idx in rep_points_gdf.index :
    n.add("Bus", f"{idx}",y=rep_points_gdf.loc[idx,'y'] , x=rep_points_gdf.loc[idx,'x'])

n.buses

carriers = ["offwind", "onwind", "solar", "coal" , "gas", "hydro", "nuclear", "oil", "hydrogen storage underground", "battery storage"]

n.madd(
    "Carrier",
    carriers,
    color=["cyan","blue","gold","grey","indianred","green","orange","chocolate","magenta", "yellowgreen"],
    co2_emissions=[costs.at[c, "CO2 intensity"] for c in carriers]
)
n.add(
    "Carrier",
    "HVAC",
    color = "brown",
    co2_emissions = 0
)

costs[['CO2 intensity']].loc['solar']

parameter
CO2 intensity    0.0
Name: solar, dtype: float64

fossils_gen = df_fossils.copy()
fossils_gen

n_fossils_gen = fossils_gen[['Coal','Gas','Nuclear']]
#n_fossils_gen = n_fossils_gen.index.astype(int)

new_index = range(1, 7)
fossils_generators = n_fossils_gen.reindex(new_index, fill_value=0)
fossils_generators = fossils_generators.rename(columns={'Gas': 'gas', 'Coal': 'coal', 'Nuclear': 'nuclear'})
fossils_generators

# Adding existing generators
generators = ["coal" , "gas", "nuclear"]

for idx in rep_points_gdf.index :
    for g in generators:
        n.add(
            "Generator",
            f"{g} {idx}",
            bus= f'{idx}',
            carrier= f"{g}",
            #capital_cost=costs.at[g, "capital_cost"],
            capital_cost=0,
            marginal_cost=costs.at[g, "marginal_cost"],
            efficiency=costs.at[g, "efficiency"],
            p_nom_extendable=False,
            p_nom = fossils_generators.loc[idx, g],
        )

n.generators

carr = "solar"
for idx in pot_max_solar.keys() : 
    n.add(
            "Generator",
            f"{carr} {idx}",
            bus=f'{idx}',
            carrier=f"{carr}",
            p_max_pu=solar_df[idx],
            capital_cost=costs.at[f"{carr}", "capital_cost"],
            marginal_cost=costs.at[f"{carr}", "marginal_cost"],
            efficiency=costs.at[f"{carr}", "efficiency"],
            p_nom_extendable=True,
            p_nom_max = pot_max_solar[idx]
        )

carr = "offwind"
for idx in pot_max_offwind.keys() : 
    n.add(
            "Generator",
            f"{carr} {idx}",
            bus=f'{idx}',
            carrier=f"{carr}",
            p_max_pu=offshore_df[idx],
            capital_cost=costs.at[f"{carr}", "capital_cost"],
            marginal_cost=costs.at[f"{carr}", "marginal_cost"],
            efficiency=costs.at[f"{carr}", "efficiency"],
            p_nom_extendable=True,
            p_nom_max = pot_max_offwind[idx]
        )
    
carr = "onwind"
for idx in pot_max_onwind.keys() : 
    n.add(
            "Generator",
            f"{carr} {idx}",
            bus=f'{idx}',
            carrier=f"{carr}",
            p_max_pu=onshore_df[idx],
            capital_cost=costs.at[f"{carr}", "capital_cost"],
            marginal_cost=costs.at[f"{carr}", "marginal_cost"],
            efficiency=costs.at[f"{carr}", "efficiency"],
            p_nom_extendable=True,
            p_nom_max = pot_max_onwind[idx]
        )

n.generators

for idx in rep_points_gdf.index:
    n.add(
        "Load",
        f"Demand {idx}",
        bus=f"{idx}",
        p_set=load.iloc[:,idx],
    )

#Adding Battery storage
for idx in rep_points_gdf.index:
    n.add(
        "StorageUnit",
        f"battery storage {idx}",
        bus=f'{idx}',
        carrier="battery storage",
        max_hours=6,
        capital_cost=costs.at["battery inverter", "capital_cost"] + 6 * costs.at["battery storage", "capital_cost"],
        efficiency_store=costs.at["battery inverter", "efficiency"],
        efficiency_dispatch=costs.at["battery inverter", "efficiency"],
        p_nom_extendable=True,
        cyclic_state_of_charge=True,
    )

#Adding Hydrogen storage
capital_costs = (
        costs.at["electrolysis", "capital_cost"] +
        costs.at["fuel cell", "capital_cost"] +
        300 * costs.at["hydrogen storage underground", "capital_cost"]
)

for idx in rep_points_gdf.index:
    n.add(
        "StorageUnit",
        f"hydrogen storage underground {idx}",
        bus=f'{idx}',
        carrier="hydrogen storage underground",
        max_hours=300,
        capital_cost=capital_costs,
        efficiency_store=costs.at["electrolysis", "efficiency"],
        efficiency_dispatch=costs.at["fuel cell", "efficiency"],
        p_nom_extendable=True,
        cyclic_state_of_charge=True,
    )

rep_points_gdf

rep_points_gdf_3035 = rep_points_gdf.to_crs(3035)
rep_points_gdf_3035

distance = rep_points_gdf_3035.loc[1].geometry.distance(rep_points_gdf_3035.loc[2].geometry)
distance

82788.24612582996

neighbours = [[1,2], [1,6], [2,3], [2,4], [2,5], [2,6], [3,5], [4,5], [4,6], [5,6]]
for pair in neighbours:
    dist = rep_points_gdf_3035.loc[pair[0]].geometry.distance(rep_points_gdf_3035.loc[pair[1]].geometry)
    n.add(
        "Link",
        f"{pair[0]}-{pair[1]}",
        bus0 = f"{pair[0]}",
        bus1 = f"{pair[1]}",
        efficiency = 1,
        marginal_cost = 0,
        p_min_pu = -1,
        length = dist * 1.5,
        capital_cost = 600 * dist *1.5,
        p_nom_extendable = True,
        carrier = "HVAC",
    )

rep_points_gdf_3035.x.iloc[0]

6.307699560864975

n.plot(margin=1, bus_sizes=0.005)

/opt/anaconda3/envs/pitie/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '

(<matplotlib.collections.PatchCollection at 0x1c60f96d0>,
 <matplotlib.collections.LineCollection at 0x191f3b850>)

n.optimize(solver_name="gurobi")

INFO:linopy.model: Solve problem using Gurobi solver

Set parameter Username

INFO:gurobipy:Set parameter Username

Academic license - for non-commercial use only - expires 2025-07-13

INFO:gurobipy:Academic license - for non-commercial use only - expires 2025-07-13
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 20/20 [00:09<00:00,  2.14it/s]
Writing continuous variables.: 100%|██████████| 8/8 [00:02<00:00,  3.87it/s]
INFO:linopy.io: Writing time: 12.21s

Read LP format model from file /private/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/linopy-problem-e_u32e8s.lp

INFO:gurobipy:Read LP format model from file /private/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/linopy-problem-e_u32e8s.lp

Reading time = 4.40 seconds

INFO:gurobipy:Reading time = 4.40 seconds

obj: 1471724 rows, 657033 columns, 2938740 nonzeros

def plot_dispatch(n, time, ymin, ymax):

    p_by_carrier = n.generators_t.p.groupby(n.generators.carrier, axis=1).sum().div(1e3)

    if not n.storage_units.empty:
        sto = n.storage_units_t.p.groupby(n.storage_units.carrier, axis=1).sum().div(1e3)
        p_by_carrier = pd.concat([p_by_carrier, sto], axis=1)

    fig, ax = plt.subplots(figsize=(25, 15))

    color = p_by_carrier.columns.map(n.carriers.color)

    p_by_carrier.where(p_by_carrier>0).loc[time].plot.area(
        ax=ax,
        linewidth=0,
        color=color,
    )
    plt.rcParams.update({'font.size': 14})

    charge = p_by_carrier.where(p_by_carrier<0).dropna(how='all', axis=1).loc[time]

    if not charge.empty:

        charge.plot.area(
            ax=ax,
            linewidth=0,
            color=charge.columns.map(n.carriers.color),
        )

    n.loads_t.p_set.sum(axis=1).loc[time].div(1e3).plot(ax=ax, c='k')
    
    plt.legend(loc=(1.01,0))
    ax.set_ylabel("GW")
    
    #set axis limit to min and max overall
    ax.set_ylim(ymin, ymax)
    # ax.set_ylim(p_by_carrier.where(p_by_carrier>0).loc[time].min(), p_by_carrier.where(p_by_carrier>0).loc[time].min().max())

plot_dispatch(n, "2013", -10, 105)

/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/ipykernel_27816/3364439340.py:3: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  p_by_carrier = n.generators_t.p.groupby(n.generators.carrier, axis=1).sum().div(1e3)
/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/ipykernel_27816/3364439340.py:6: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  sto = n.storage_units_t.p.groupby(n.storage_units.carrier, axis=1).sum().div(1e3)
/opt/anaconda3/envs/pitie/lib/python3.11/site-packages/pandas/plotting/_matplotlib/core.py:1794: UserWarning: Attempting to set identical low and high ylims makes transformation singular; automatically expanding.
  ax.set_ylim(None, 0)

n1 = n.copy()

n1.add(
    "GlobalConstraint",
    "CO2Limit",
    carrier_attribute="co2_emissions",
    sense="<=",
    constant=0,
)

# for idx in rep_points_gdf.index :
#     n1.remove("Generator",f"gas {idx}")
#     n1.remove("Generator",f"coal {idx}")

n1.generators

n1.optimize(solver_name="gurobi")

INFO:linopy.model: Solve problem using Gurobi solver

Set parameter Username

INFO:gurobipy:Set parameter Username

Academic license - for non-commercial use only - expires 2025-07-13

INFO:gurobipy:Academic license - for non-commercial use only - expires 2025-07-13
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 21/21 [00:20<00:00,  1.03it/s]
Writing continuous variables.: 100%|██████████| 8/8 [00:03<00:00,  2.17it/s]
INFO:linopy.io: Writing time: 24.72s

Read LP format model from file /private/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/linopy-problem-407vl4sn.lp

INFO:gurobipy:Read LP format model from file /private/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/linopy-problem-407vl4sn.lp

Reading time = 4.55 seconds

INFO:gurobipy:Reading time = 4.55 seconds

obj: 1471725 rows, 657033 columns, 3043860 nonzeros

plot_dispatch(n1, "2013", -10, 105)

/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/ipykernel_27816/3364439340.py:3: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  p_by_carrier = n.generators_t.p.groupby(n.generators.carrier, axis=1).sum().div(1e3)
/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/ipykernel_27816/3364439340.py:6: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  sto = n.storage_units_t.p.groupby(n.storage_units.carrier, axis=1).sum().div(1e3)
/opt/anaconda3/envs/pitie/lib/python3.11/site-packages/pandas/plotting/_matplotlib/core.py:1794: UserWarning: Attempting to set identical low and high ylims makes transformation singular; automatically expanding.
  ax.set_ylim(None, 0)

buses_marginal_price_0 = n.buses_t.marginal_price.copy()
buses_marginal_price_1 = n1.buses_t.marginal_price.copy()

buses_marginal_price_0

buses_marginal_price_1

n1.generators.p_nom_opt

Generator
coal 1        1600.000000
gas 1         4306.000000
nuclear 1        0.000000
coal 2         590.000000
gas 2          349.000000
nuclear 2        0.000000
coal 3           0.000000
gas 3            0.000000
nuclear 3        0.000000
coal 4         685.000000
gas 4         1025.000000
nuclear 4        0.000000
coal 5        2140.000000
gas 5         1753.000000
nuclear 5      485.000000
coal 6           0.000000
gas 6         1879.000000
nuclear 6        0.000000
solar 1       4789.590677
solar 2      18207.071100
solar 3      42201.000000
solar 4      60556.947182
solar 5      38138.223007
offwind 6        0.000000
onwind 1      5924.845681
onwind 2     10357.650000
onwind 3      3581.100000
onwind 4      2205.450000
onwind 5      3563.100000
Name: p_nom_opt, dtype: float64

n.storage_units.p_nom_opt

StorageUnit
battery storage 1                    0.000000
battery storage 2                 1896.206441
battery storage 3                 9995.955392
battery storage 4                 5180.875353
battery storage 5                    0.000000
battery storage 6                    0.000000
hydrogen storage underground 1       0.000000
hydrogen storage underground 2    1474.031973
hydrogen storage underground 3    3462.055581
hydrogen storage underground 4     873.125349
hydrogen storage underground 5       0.000000
hydrogen storage underground 6       0.000000
Name: p_nom_opt, dtype: float64

n.storage_units.max_hours

StorageUnit
battery storage 1                   6.0
battery storage 2                   6.0
battery storage 3                   6.0
battery storage 4                   6.0
battery storage 5                   6.0
battery storage 6                   6.0
hydrogen storage underground 1    300.0
hydrogen storage underground 2    300.0
hydrogen storage underground 3    300.0
hydrogen storage underground 4    300.0
hydrogen storage underground 5    300.0
hydrogen storage underground 6    300.0
Name: max_hours, dtype: float64

n.storage_units.p_nom_opt.div(1e3) * n.storage_units.max_hours  # GWh

StorageUnit
battery storage 1                    0.000000
battery storage 2                   11.377239
battery storage 3                   59.975732
battery storage 4                   31.085252
battery storage 5                    0.000000
battery storage 6                    0.000000
hydrogen storage underground 1       0.000000
hydrogen storage underground 2     442.209592
hydrogen storage underground 3    1038.616674
hydrogen storage underground 4     261.937605
hydrogen storage underground 5       0.000000
hydrogen storage underground 6       0.000000
dtype: float64

def system_cost(n):
    tsc = n.statistics.capex() + n.statistics.opex()
    gdf_combined = n.statistics.capex().add(n.statistics.opex(), fill_value=0)
    return gdf_combined.droplevel(0)/1e6  # million €/a

costs0 = system_cost(n)
costs1 = system_cost(n1)
#costs = costs.fillna(0)
#costs0 = costs0.dropna()
#costs1 = costs1.dropna()
#costs0

#system_cost(n1).plot.pie(figsize=(2, 2))

type(costs0)

pandas.core.series.Series

def func(pct, allvals):
    return "{:.2f}%".format(pct)

# Create the pie chart
fig, ax = plt.subplots(figsize=(6, 6))
wedges, texts, autotexts = ax.pie(costs0, labels=costs0.index, autopct=lambda pct: func(pct, costs0), startangle=140, textprops=dict(color="w"))

# Format the text on the pie chart
for text in autotexts:
    text.set_color('black')
    text.set_weight('bold')

# Create custom legend with absolute values
legend_labels = [f"{label}: €{value:.1f}M" for label, value in zip(costs0.index, costs0)]
ax.legend(wedges, legend_labels, title="Cost Categories", loc="center left", bbox_to_anchor=(1, 0, 0.5, 1))

ax.set_title("System Costs Distribution (Capx + Opex) for Emissions Scenario")

total_cost = costs0.sum()
plt.text(0, -1.2, f"Total Cost: €{total_cost:.1f}M", ha='center', va='center', fontsize=12, fontweight='bold')

plt.show()

fig, ax = plt.subplots(figsize=(6, 6))
wedges, texts, autotexts = ax.pie(costs1, labels=costs1.index, autopct=lambda pct: func(pct, costs1), startangle=140, textprops=dict(color="w"))

# Format the text on the pie chart
for text in autotexts:
    text.set_color('black')
    text.set_weight('bold')

# Create custom legend with absolute values
legend_labels = [f"{label}: €{value:.2f}M" for label, value in zip(costs1.index, costs1)]
ax.legend(wedges, legend_labels, title="Cost Categories", loc="center left", bbox_to_anchor=(1, 0, 0.5, 1))

total_cost = costs1.sum()
plt.text(0, -1.2, f"Total Cost: €{total_cost:.1f}M", ha='center', va='center', fontsize=12, fontweight='bold')

ax.set_title("System Costs Distribution (Capx + Opex) for No Emissions Scenario")
plt.show()

print("costs scenario Emissions M€", costs0.sum())
print("costs scenario No Emissions M€", costs1.sum())

costs scenario Emissions M€ 13165.878462525432
costs scenario No Emissions M€ 24235.495203921328

n.statistics.capex()

component    carrier                     
Generator    onwind                          1.670786e+09
             solar                           4.791859e+09
StorageUnit  battery storage                 2.377848e+09
             hydrogen storage underground    3.057820e+09
dtype: float64

n.statistics.opex().droplevel(0)

carrier
coal       3.776760e+08
gas        7.908262e+08
nuclear    5.953678e+07
onwind     3.888908e+07
solar      6.378157e+05
dtype: float64

cap1 = n1.generators.groupby("carrier").p_nom_opt.sum()
cap1 = cap1/1e3

def func(pct, allvals):
    return "{:.1f}%".format(pct)

# Create the pie chart
fig, ax = plt.subplots(figsize=(6, 6))
wedges, texts, autotexts = ax.pie(cap1, labels=cap1.index, autopct=lambda pct: func(pct, cap1), startangle=140, textprops=dict(color="w"))

# Format the text on the pie chart
for text in autotexts:
    text.set_color('black')
    text.set_weight('bold')

# Create custom legend with absolute values
legend_labels = [f"{label}: {value:.1f}GW" for label, value in zip(cap1.index, cap1)]
ax.legend(wedges, legend_labels, title="Carrier", loc="center left", bbox_to_anchor=(1, 0, 0.5, 1))

ax.set_title("System Capacities for No CO2 Emissions Scenario")
plt.show()

cap0 = n.generators.groupby("carrier").p_nom_opt.sum()
cap0 = cap0/1e3

# Create the pie chart
fig, ax = plt.subplots(figsize=(6, 6))
wedges, texts, autotexts = ax.pie(cap0, labels=cap0.index, autopct=lambda pct: func(pct, cap0), startangle=140, textprops=dict(color="w"))

# Format the text on the pie chart
for text in autotexts:
    text.set_color('black')
    text.set_weight('bold')

# Create custom legend with absolute values
legend_labels = [f"{label}: {value:.1f}GW" for label, value in zip(cap0.index, cap0)]
ax.legend(wedges, legend_labels, title="Carrier", loc="center left", bbox_to_anchor=(1, 0, 0.5, 1))

ax.set_title("System Capacities for CO2 Emissions Scenario")
plt.show()

test = n.storage_units_t.state_of_charge.groupby([n.storage_units.bus], axis=1).sum()
regioncolors = ["orange", "skyblue", "brown", "grey", "lightgreen", "purple"]
region_names = test.columns  # Assuming region_names should match the columns in test
i = 0
# test = n.storage_units_t.state_of_charge.groupby([n.storage_units.bus], axis=1).sum()

fig, ax = plt.subplots(figsize=(15, 10))

for region in region_names:
    if region in test.columns:
        x = test[region] / 1e3
        x.plot(color=regioncolors[i % len(regioncolors)], label=region)
        i += 1
    else:
        print(f"Region {region} not found in test columns.")

plt.title('Storage State of Charge for each Region with CO2 Emissions')
plt.legend(loc='upper left')
plt.ylabel('State of Charge (GWh)')
plt.xlabel('2013')
plt.show()


# The state of charge (SOC) of a cell denotes the capacity that is currently available as a function of the rated capacity. The value of the SOC varies between 0% and 100%.

/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/ipykernel_27816/2483751166.py:1: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  test = n.storage_units_t.state_of_charge.groupby([n.storage_units.bus], axis=1).sum()

test = n1.storage_units_t.state_of_charge.groupby([n1.storage_units.bus], axis=1).sum()
regioncolors = ["orange", "skyblue", "brown", "grey", "lightgreen", "purple"]
region_names = test.columns  # Assuming region_names should match the columns in test
i = 0
# test = n.storage_units_t.state_of_charge.groupby([n.storage_units.bus], axis=1).sum()

fig, ax = plt.subplots(figsize=(15, 10))

for region in region_names:
    if region in test.columns:
        x = test[region] / 1e3
        x.plot(color=regioncolors[i % len(regioncolors)], label=region)
        i += 1
    else:
        print(f"Region {region} not found in test columns.")

plt.title('Storage State of Charge for each Region without CO2 Emissions')
plt.legend(loc='upper left')
plt.ylabel('State of Charge (GWh)')
plt.xlabel('2013')
plt.show()

/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/ipykernel_27816/1214001169.py:1: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  test = n1.storage_units_t.state_of_charge.groupby([n1.storage_units.bus], axis=1).sum()

price_duration=n.buses_t.marginal_price
price_duration=price_duration.mean(axis=1)
price_duration=price_duration.sort_values(ascending=False).reset_index(drop=True)
price_duration.plot(ylabel='Cost [Eur]', xlabel='frequency [1/year]', title='Price Duration Curve Emissions Scenario')
plt.ylim(0,2000)

(0.0, 2000.0)

price_duration=n1.buses_t.marginal_price
price_duration=price_duration.mean(axis=1)
price_duration=price_duration.sort_values(ascending=False).reset_index(drop=True)
price_duration.plot(ylabel='Cost [Eur]', xlabel='frequency [1/year]]', title='Price Duration Curve No Emissions Scenario')
plt.ylim(0,2000)

(0.0, 2000.0)

avg_price_r = n.buses_t.marginal_price.mean()

colors = ['blue', 'green', 'red', 'purple', 'orange', 'cyan']

# Plotting the bar graph
fig, ax = plt.subplots(figsize=(10, 5))
avg_price_r.plot(kind='bar', color=colors, ax=ax)

# Adding title and labels
ax.set_title('Average Price by Region Emissions Scenario')
ax.set_xlabel('Region')
ax.set_ylabel('Average Price [€/MWh]')

# Adding legend
handles = [plt.Rectangle((0,0),1,1, color=color) for color in colors]
labels = avg_price_r.index
ax.legend(handles, labels, title='Regions', loc='best')

# Display the plot
plt.show()

avg_price_r1 = n1.buses_t.marginal_price.mean()

colors = ['blue', 'green', 'red', 'purple', 'orange', 'cyan']

# Plotting the bar graph
fig, ax = plt.subplots(figsize=(10, 5))
avg_price_r1.plot(kind='bar', color=colors, ax=ax)

# Adding title and labels
ax.set_title('Average Price by Region No Emissions Scenario')
ax.set_xlabel('Region')
ax.set_ylabel('Average Price [€/MWh]')

# Adding legend
handles = [plt.Rectangle((0,0),1,1, color=color) for color in colors]
labels = avg_price_r.index
ax.legend(handles, labels, title='Regions', loc='best')

# Display the plot
plt.show()

fig = plt.figure(figsize=(7,7))
ax = plt.axes(projection=ccrs.PlateCarree())

norm = plt.Normalize(vmin=n.buses_t.marginal_price.mean().min(), vmax=n.buses_t.marginal_price.mean().max()) # €/MWh

n.plot(
    margin=0.7,
    ax=ax,
    bus_colors=n.buses_t.marginal_price.mean(),
    bus_cmap='plasma',
    bus_norm=norm,
)

sm = plt.cm.ScalarMappable(cmap='plasma', norm=norm)
sm.set_array([])

ax.set_title("Average Marginal Price Distribution Emissions Scenario")

cbar = plt.colorbar(sm, ax=ax, label="Average Marginal Price [€/MWh]", shrink=0.6)

/opt/anaconda3/envs/pitie/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '

fig = plt.figure(figsize=(7,7))
ax = plt.axes(projection=ccrs.PlateCarree())

norm = plt.Normalize(vmin=n1.buses_t.marginal_price.mean().min(), vmax=n1.buses_t.marginal_price.mean().max()) # €/MWh

n1.plot(
    margin=0.7,
    ax=ax,
    bus_colors=n1.buses_t.marginal_price.mean(),
    bus_cmap='plasma',
    bus_norm=norm,
)

sm = plt.cm.ScalarMappable(cmap='plasma', norm=norm)
sm.set_array([])

ax.set_title("Average Marginal Price Distribution No Emissions Scenario")

cbar = plt.colorbar(sm, ax=ax, label="Average Marginal Price [€/MWh]", shrink=0.6)

/opt/anaconda3/envs/pitie/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '

# CHECKING THE GENERATION OF THE POWERPLANTS AFTER OPTIMIZATION 
solar_generators = n1.generators[n1.generators.carrier == 'solar']

# Obtenez les générations optimisées pour ces générateurs solaires
solar_generation = n1.generators_t.p[solar_generators.index]

solar_generation_df = pd.DataFrame(solar_generation)

# Plot the generation of solar power plants
solar_generation_df.plot(figsize=(10, 6))
plt.title('Generation of Solar Power Plants After Optimization')
plt.xlabel('Time')
plt.ylabel('Generation (MW)')
plt.legend(title='Solar Power Plants')
plt.show()

shadow_price_CO2 = n1.global_constraints.loc['CO2Limit', 'mu']

# I am surprised by the value, is the shadow price the amount of money it would cost to reduce the CO2 emissions by one unit? 
# Given that we set up the constraint at 0, does it make sense ? 

print("In the Zero Emissions Scenario, the Shadow Price of CO2 is: ", f"{np.abs(shadow_price_CO2):.2f}", "€/tCO2")

In the Zero Emissions Scenario, the Shadow Price of CO2 is:  32736.40 €/tCO2

capacity_factors = n.generators_t.p_max_pu

# Obtenez la génération effective pour chaque générateur
generation_effective = n.generators_t.p

# Filtrer les générateurs par type
solar_generators = n.generators[n.generators.carrier == 'solar']
onwind_generators = n.generators[n.generators.carrier == 'onwind']
offwind_generators = n.generators[n.generators.carrier == 'offwind']

# Utiliser les capacités optimisées p_nom_opt
solar_p_nom_opt = n.generators.loc[solar_generators.index, 'p_nom_opt']
onwind_p_nom_opt = n.generators.loc[onwind_generators.index, 'p_nom_opt']
offwind_p_nom_opt = n.generators.loc[offwind_generators.index, 'p_nom_opt']

# Calculer le curtailment pour chaque type de générateur à chaque instant de temps
curtailment_solar = capacity_factors[solar_generators.index]*solar_p_nom_opt - generation_effective[solar_generators.index]
curtailment_onwind = capacity_factors[onwind_generators.index]*onwind_p_nom_opt - generation_effective[onwind_generators.index]
curtailment_offwind = capacity_factors[offwind_generators.index]*offwind_p_nom_opt - generation_effective[offwind_generators.index]

fig, axs = plt.subplots(3, 1, figsize=(12, 18))
fig.suptitle('Curtailments for the Emissions Scenario', fontsize=16)

curtailment_solar.sum(axis=1).plot(ax=axs[0], label='Solar Curtailment', color='orange')
axs[0].set_title('Curtailment of Solar Generators')
axs[0].set_xlabel('Time')
axs[0].set_ylabel('Curtailment (MW)')
axs[0].legend()

curtailment_onwind.sum(axis=1).plot(ax=axs[1], label='Onshore Wind Curtailment', color='blue')
axs[1].set_title('Curtailment of Onshore Wind Generators')
axs[1].set_xlabel('Time')
axs[1].set_ylabel('Curtailment (MW)')
axs[1].legend()

curtailment_offwind.sum(axis=1).plot(ax=axs[2], label='Offshore Wind Curtailment', color='green')
axs[2].set_title('Curtailment of Offshore Wind Generators')
axs[2].set_xlabel('Time')
axs[2].set_ylabel('Curtailment (MW)')
axs[2].legend()

plt.tight_layout()
plt.show()

capacity_factors = n1.generators_t.p_max_pu

# Obtenez la génération effective pour chaque générateur
generation_effective = n1.generators_t.p

# Filtrer les générateurs par type
solar_generators = n1.generators[n1.generators.carrier == 'solar']
onwind_generators = n1.generators[n1.generators.carrier == 'onwind']
offwind_generators = n1.generators[n1.generators.carrier == 'offwind']

# Utiliser les capacités optimisées p_nom_opt
solar_p_nom_opt = n1.generators.loc[solar_generators.index, 'p_nom_opt']
onwind_p_nom_opt = n1.generators.loc[onwind_generators.index, 'p_nom_opt']
offwind_p_nom_opt = n1.generators.loc[offwind_generators.index, 'p_nom_opt']

# Calculer le curtailment pour chaque type de générateur à chaque instant de temps
curtailment_solar = capacity_factors[solar_generators.index]*solar_p_nom_opt - generation_effective[solar_generators.index]
curtailment_onwind = capacity_factors[onwind_generators.index]*onwind_p_nom_opt - generation_effective[onwind_generators.index]
curtailment_offwind = capacity_factors[offwind_generators.index]*offwind_p_nom_opt - generation_effective[offwind_generators.index]

fig, axs = plt.subplots(3, 1, figsize=(12, 18))
fig.suptitle('Curtailments for the No Emissions Scenario', fontsize=16)

curtailment_solar.sum(axis=1).plot(ax=axs[0], label='Solar Curtailment', color='orange')
axs[0].set_title('Curtailment of Solar Generators')
axs[0].set_xlabel('Time')
axs[0].set_ylabel('Curtailment (MW)')
axs[0].legend()

curtailment_onwind.sum(axis=1).plot(ax=axs[1], label='Onshore Wind Curtailment', color='blue')
axs[1].set_title('Curtailment of Onshore Wind Generators')
axs[1].set_xlabel('Time')
axs[1].set_ylabel('Curtailment (MW)')
axs[1].legend()

curtailment_offwind.sum(axis=1).plot(ax=axs[2], label='Offshore Wind Curtailment', color='green')
axs[2].set_title('Curtailment of Offshore Wind Generators')
axs[2].set_xlabel('Time')
axs[2].set_ylabel('Curtailment (MW)')
axs[2].legend()

plt.tight_layout()
plt.show()

n2 = n1.copy()

for idx in rep_points_gdf.index :
    n2.remove("Generator",f"nuclear {idx}")

n2.add(
            "Generator",
            f"nuclear 5",
            bus= "5",
            carrier= "nuclear",
            capital_cost=0,
            marginal_cost=costs.at["nuclear", "marginal_cost"],
            efficiency=costs.at[g, "efficiency"],
            p_nom_extendable=False,
            p_nom = fossils_generators.loc[5, "nuclear"],
        )

n2.generators

n2.carriers['co2_emissions']

Carrier
offwind                         0.0000
onwind                          0.0000
solar                           0.0000
coal                            0.3361
gas                             0.1980
hydro                           0.0000
nuclear                         0.0000
oil                             0.2571
hydrogen storage underground    0.0000
battery storage                 0.0000
HVAC                            0.0000
Name: co2_emissions, dtype: float64

n2.generators['capital_cost']

Generator
coal 1            0.000000
gas 1             0.000000
coal 2            0.000000
gas 2             0.000000
coal 3            0.000000
gas 3             0.000000
coal 4            0.000000
gas 4             0.000000
coal 5            0.000000
gas 5             0.000000
coal 6            0.000000
gas 6             0.000000
solar 1       63114.923985
solar 2       63114.923985
solar 3       63114.923985
solar 4       63114.923985
solar 5       63114.923985
offwind 6    184567.393550
onwind 1     107441.457021
onwind 2     107441.457021
onwind 3     107441.457021
onwind 4     107441.457021
onwind 5     107441.457021
nuclear 5         0.000000
Name: capital_cost, dtype: float64

costs.loc["nuclear"]

parameter
C in fuel                                 NaN
C stored                                  NaN
CO2 intensity                    0.000000e+00
CO2 stored                                NaN
FOM                              1.270000e+00
Motor size                                NaN
VOM                              3.546400e+00
ammonia-input                             NaN
c_b                                       NaN
c_v                                       NaN
capacity                                  NaN
capture rate                              NaN
capture_rate                              NaN
carbondioxide-input                       NaN
carbondioxide-output                      NaN
commodity                                 NaN
compression-electricity-input             NaN
compression-heat-output                   NaN
discount rate                    7.000000e-02
district heat surcharge                   NaN
district heat-input                       NaN
efficiency                       3.260000e-01
efficiency-biomass                        NaN
efficiency-electricity                    NaN
efficiency-heat                           NaN
efficiency-hydrogen                       NaN
efficiency-tot                            NaN
electricity-input                         NaN
fuel                             3.412200e+00
hbi-input                                 NaN
heat-input                                NaN
heat-output                               NaN
hydrogen-input                            NaN
investment                       8.594135e+06
lifetime                         4.000000e+01
lohc-input                                NaN
methane-input                             NaN
methanol-input                            NaN
min_fill_level                            NaN
naphtha-input                             NaN
nitrogen-input                            NaN
ore-input                                 NaN
p_nom_ratio                               NaN
pelletizing cost                          NaN
marginal_cost                    1.401327e+01
capital_cost                     7.537842e+05
Name: nuclear, dtype: float64

costs[['capital_cost']].loc['nuclear']

parameter
capital_cost    753784.215297
Name: nuclear, dtype: float64

costs_df.index.get_level_values(0)

Index(['Ammonia cracker', 'Ammonia cracker', 'Ammonia cracker',
       'Ammonia cracker', 'BEV Bus city', 'BEV Bus city', 'BEV Bus city',
       'BEV Bus city', 'BEV Bus city', 'BEV Bus city',
       ...
       'waste CHP CC', 'waste CHP CC', 'waste CHP CC', 'waste CHP CC',
       'waste CHP CC', 'waste CHP CC', 'waste CHP CC', 'waste CHP CC',
       'water tank charger', 'water tank discharger'],
      dtype='object', name='technology', length=1033)

gas_costs = costs_df.xs('gas', level=0)
gas_costs

costs_df[costs_df.index.get_level_values(0) == 'nuclear']

values_to_keep = ["nuclear", "offwind", "onwind", "solar"]

# Filtrer le DataFrame pour ne conserver que les lignes correspondant à ces valeurs
costs[costs.index.isin(values_to_keep)]

inf, sup = 2.5e6, 10e6

values = np.linspace(inf, sup, 4)
nuclear_cc = values.tolist()
nuclear_cc

[2500000.0, 5000000.0, 7500000.0, 10000000.0]

#p_nom = fossils_generators.loc[idx, g]
fossils_generators

carr = "nuclear"
for idx in rep_points_gdf.index : 
    n2.add(
            "Generator",
            f"{carr} add {idx}",
            bus=f'{idx}',
            carrier=f"{carr}",
            capital_cost= nuclear_cc[0],
            marginal_cost=costs.at[f"{carr}", "marginal_cost"],
            efficiency=costs.at[f"{carr}", "efficiency"],
            p_nom_extendable=True,
            p_nom_max = fossils_generators.loc[5, "nuclear"]
        )

n2.generators['capital_cost']

Generator
coal 1           0.000000e+00
gas 1            0.000000e+00
coal 2           0.000000e+00
gas 2            0.000000e+00
coal 3           0.000000e+00
gas 3            0.000000e+00
coal 4           0.000000e+00
gas 4            0.000000e+00
coal 5           0.000000e+00
gas 5            0.000000e+00
coal 6           0.000000e+00
gas 6            0.000000e+00
solar 1          6.311492e+04
solar 2          6.311492e+04
solar 3          6.311492e+04
solar 4          6.311492e+04
solar 5          6.311492e+04
offwind 6        1.845674e+05
onwind 1         1.074415e+05
onwind 2         1.074415e+05
onwind 3         1.074415e+05
onwind 4         1.074415e+05
onwind 5         1.074415e+05
nuclear 5        0.000000e+00
nuclear add 1    2.500000e+06
nuclear add 2    2.500000e+06
nuclear add 3    2.500000e+06
nuclear add 4    2.500000e+06
nuclear add 5    2.500000e+06
nuclear add 6    2.500000e+06
Name: capital_cost, dtype: float64

n2.optimize(solver_name="gurobi")

INFO:linopy.model: Solve problem using Gurobi solver

Set parameter Username

INFO:gurobipy:Set parameter Username

Academic license - for non-commercial use only - expires 2025-07-13

INFO:gurobipy:Academic license - for non-commercial use only - expires 2025-07-13
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 21/21 [00:14<00:00,  1.47it/s]
Writing continuous variables.: 100%|██████████| 8/8 [00:02<00:00,  3.70it/s]
INFO:linopy.io: Writing time: 17.3s

Read LP format model from file /private/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/linopy-problem-_u5i728e.lp

INFO:gurobipy:Read LP format model from file /private/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/linopy-problem-_u5i728e.lp

Reading time = 5.18 seconds

INFO:gurobipy:Reading time = 5.18 seconds

obj: 1489257 rows, 665799 columns, 3122712 nonzeros

n3 = n2.copy()

ncc = nuclear_cc[1] 

for x in range(1, 7):
    generator_name = f"nuclear add {x}"
    if generator_name in n3.generators.index:
        n3.generators.at[generator_name, 'capital_cost'] = ncc

n3.generators

n3.optimize(solver_name="gurobi")

INFO:linopy.model: Solve problem using Gurobi solver

Set parameter Username

INFO:gurobipy:Set parameter Username

Academic license - for non-commercial use only - expires 2025-07-13

INFO:gurobipy:Academic license - for non-commercial use only - expires 2025-07-13
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 21/21 [00:14<00:00,  1.42it/s]
Writing continuous variables.: 100%|██████████| 8/8 [00:02<00:00,  3.84it/s]
INFO:linopy.io: Writing time: 17.87s

Read LP format model from file /private/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/linopy-problem-yfkvjzkb.lp

INFO:gurobipy:Read LP format model from file /private/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/linopy-problem-yfkvjzkb.lp

Reading time = 5.61 seconds

INFO:gurobipy:Reading time = 5.61 seconds

obj: 1489257 rows, 665799 columns, 3122712 nonzeros

n4 = n2.copy()

ncc = nuclear_cc[2] 

for x in range(1, 7):
    generator_name = f"nuclear add {x}"
    if generator_name in n4.generators.index:
        n4.generators.at[generator_name, 'capital_cost'] = ncc

n4.generators['capital_cost']

Generator
coal 1           0.000000e+00
gas 1            0.000000e+00
coal 2           0.000000e+00
gas 2            0.000000e+00
coal 3           0.000000e+00
gas 3            0.000000e+00
coal 4           0.000000e+00
gas 4            0.000000e+00
coal 5           0.000000e+00
gas 5            0.000000e+00
coal 6           0.000000e+00
gas 6            0.000000e+00
solar 1          6.311492e+04
solar 2          6.311492e+04
solar 3          6.311492e+04
solar 4          6.311492e+04
solar 5          6.311492e+04
offwind 6        1.845674e+05
onwind 1         1.074415e+05
onwind 2         1.074415e+05
onwind 3         1.074415e+05
onwind 4         1.074415e+05
onwind 5         1.074415e+05
nuclear 5        0.000000e+00
nuclear add 1    7.500000e+06
nuclear add 2    7.500000e+06
nuclear add 3    7.500000e+06
nuclear add 4    7.500000e+06
nuclear add 5    7.500000e+06
nuclear add 6    7.500000e+06
Name: capital_cost, dtype: float64

n4.optimize(solver_name="gurobi")

INFO:linopy.model: Solve problem using Gurobi solver

Set parameter Username

INFO:gurobipy:Set parameter Username

Academic license - for non-commercial use only - expires 2025-07-13

INFO:gurobipy:Academic license - for non-commercial use only - expires 2025-07-13
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 21/21 [00:13<00:00,  1.62it/s]
Writing continuous variables.: 100%|██████████| 8/8 [00:01<00:00,  4.07it/s]
INFO:linopy.io: Writing time: 15.84s

Read LP format model from file /private/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/linopy-problem-zrjdxt3o.lp

INFO:gurobipy:Read LP format model from file /private/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/linopy-problem-zrjdxt3o.lp

Reading time = 5.17 seconds

INFO:gurobipy:Reading time = 5.17 seconds

obj: 1489257 rows, 665799 columns, 3122712 nonzeros

n5 = n2.copy()

ncc = nuclear_cc[3] 

for x in range(1, 7):
    generator_name = f"nuclear add {x}"
    if generator_name in n5.generators.index:
        n5.generators.at[generator_name, 'capital_cost'] = ncc

n5.generators['capital_cost']

Generator
coal 1           0.000000e+00
gas 1            0.000000e+00
coal 2           0.000000e+00
gas 2            0.000000e+00
coal 3           0.000000e+00
gas 3            0.000000e+00
coal 4           0.000000e+00
gas 4            0.000000e+00
coal 5           0.000000e+00
gas 5            0.000000e+00
coal 6           0.000000e+00
gas 6            0.000000e+00
solar 1          6.311492e+04
solar 2          6.311492e+04
solar 3          6.311492e+04
solar 4          6.311492e+04
solar 5          6.311492e+04
offwind 6        1.845674e+05
onwind 1         1.074415e+05
onwind 2         1.074415e+05
onwind 3         1.074415e+05
onwind 4         1.074415e+05
onwind 5         1.074415e+05
nuclear 5        0.000000e+00
nuclear add 1    1.000000e+07
nuclear add 2    1.000000e+07
nuclear add 3    1.000000e+07
nuclear add 4    1.000000e+07
nuclear add 5    1.000000e+07
nuclear add 6    1.000000e+07
Name: capital_cost, dtype: float64

n5.optimize(solver_name = "gurobi")

INFO:linopy.model: Solve problem using Gurobi solver

Set parameter Username

INFO:gurobipy:Set parameter Username

Academic license - for non-commercial use only - expires 2025-07-13

INFO:gurobipy:Academic license - for non-commercial use only - expires 2025-07-13
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 21/21 [00:23<00:00,  1.14s/it]
Writing continuous variables.: 100%|██████████| 8/8 [00:02<00:00,  3.61it/s]
INFO:linopy.io: Writing time: 27.58s

Read LP format model from file /private/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/linopy-problem-9bvvpxub.lp

INFO:gurobipy:Read LP format model from file /private/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/linopy-problem-9bvvpxub.lp

Reading time = 5.54 seconds

INFO:gurobipy:Reading time = 5.54 seconds

obj: 1489257 rows, 665799 columns, 3122712 nonzeros

def get_generation_by_carrier(network):
    generation = network.generators_t.p.groupby(network.generators.carrier, axis=1).sum().sum()
    return generation

# Get generation by carrier for each network
gen_n2 = get_generation_by_carrier(n2)
gen_n3 = get_generation_by_carrier(n3)
gen_n4 = get_generation_by_carrier(n4)
gen_n5 = get_generation_by_carrier(n5)

# Define colors for the pie charts
colors = ['#ff9999','#66b3ff','#99ff99','#ffcc99','#c2c2f0','#ffb3e6']

# Create the figure and subplots for the pie charts
fig, axs = plt.subplots(2, 2, figsize=(14, 10))

# Main title for the figure
fig.suptitle('Comparison of Generation by Carrier After Optimization', fontsize=16)

# Function to plot a pie chart
def plot_pie(ax, generation, title):
    wedges, texts = ax.pie(generation, labels=None, startangle=90, colors=colors)
    ax.set_title(title)
    ax.axis('equal')  # Ensure the pie chart is circular
    # Create legend labels with values and percentages
    total_gen = generation.sum()
    legend_labels = [f'{carrier} ({gen/1e3:.2f} GW, {gen/total_gen:.1%})' for carrier, gen in generation.items()]
    ax.legend(wedges, legend_labels, title="Carriers", loc="center left", bbox_to_anchor=(1, 0, 0.5, 1))

# Plot the pie charts for each network
plot_pie(axs[0, 0], gen_n2, 'Nuclear Capital Cost : 2500€/kw')
plot_pie(axs[0, 1], gen_n3, 'Nuclear Capital Cost : 5000€/kw')
plot_pie(axs[1, 0], gen_n4, 'Nuclear Capital Cost : 7500€/kw')
plot_pie(axs[1, 1], gen_n5, 'Nuclear Capital Cost : 10000€/kw')

# Adjust the layout to make room for the legend
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()

/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/ipykernel_27816/2256937941.py:2: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  generation = network.generators_t.p.groupby(network.generators.carrier, axis=1).sum().sum()
/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/ipykernel_27816/2256937941.py:2: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  generation = network.generators_t.p.groupby(network.generators.carrier, axis=1).sum().sum()
/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/ipykernel_27816/2256937941.py:2: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  generation = network.generators_t.p.groupby(network.generators.carrier, axis=1).sum().sum()
/var/folders/fw/f2yz5zns65jgfw709v299lb00000gp/T/ipykernel_27816/2256937941.py:2: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  generation = network.generators_t.p.groupby(network.generators.carrier, axis=1).sum().sum()

region	1	2	3	4	5
time
2013-01-01 00:00:00	0.998975	0.997925	0.998468	0.998050	0.998297
2013-01-01 01:00:00	0.981815	0.995112	0.997072	0.988215	0.993124
2013-01-01 02:00:00	0.952146	0.981075	0.988307	0.946122	0.970967
2013-01-01 03:00:00	0.886245	0.916817	0.931022	0.767314	0.861064
2013-01-01 04:00:00	0.655116	0.762757	0.822314	0.523018	0.680503
...	...	...	...	...	...
2013-12-31 19:00:00	0.944623	0.898349	0.790733	0.951875	0.942194
2013-12-31 20:00:00	0.940086	0.876553	0.778794	0.944519	0.938032
2013-12-31 21:00:00	0.920623	0.837225	0.773567	0.938595	0.911153
2013-12-31 22:00:00	0.751153	0.649190	0.676206	0.857340	0.764536
2013-12-31 23:00:00	0.766319	0.694857	0.699070	0.887988	0.692507

	country	country_long	name	gppd_idnr	capacity_mw	latitude	longitude	primary_fuel	other_fuel1	other_fuel2	...	estimated_generation_gwh_2014	estimated_generation_gwh_2015	estimated_generation_gwh_2016	estimated_generation_gwh_2017	estimated_generation_note_2013	estimated_generation_note_2014	estimated_generation_note_2015	estimated_generation_note_2016	estimated_generation_note_2017	geometry
17685	NLD	Netherlands	Ameland	WKS0066277	6.0	53.4503	5.6734	Solar	NaN	NaN	...	8.95	8.98	9.00	8.99	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	POINT (5.67340 53.45030)
17686	NLD	Netherlands	Andijk	WKS0072735	15.0	52.7370	5.2290	Solar	NaN	NaN	...	22.77	22.51	22.53	22.52	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	POINT (5.22900 52.73700)
17687	NLD	Netherlands	Beabuorren	WRI1005281	10.4	53.0417	5.5014	Wind	NaN	NaN	...	22.29	23.53	24.76	22.35	WIND-V1	WIND-V1	WIND-V1	WIND-V1	WIND-V1	POINT (5.50140 53.04170)
17688	NLD	Netherlands	Borssele Nuclear Power Station	WRI1019300	485.0	51.4312	3.7174	Nuclear	NaN	NaN	...	NaN	NaN	NaN	3203.82	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	CAPACITY-FACTOR-V1	POINT (3.71740 51.43120)
17689	NLD	Netherlands	Budel-Dorplein	WKS0071714	43.9	51.2360	5.6150	Solar	NaN	NaN	...	62.58	62.52	63.02	62.02	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	POINT (5.61500 51.23600)

	country	country_long	name	gppd_idnr	capacity_mw	latitude	longitude	primary_fuel	other_fuel1	other_fuel2	...	estimated_generation_gwh_2014	estimated_generation_gwh_2015	estimated_generation_gwh_2016	estimated_generation_gwh_2017	estimated_generation_note_2013	estimated_generation_note_2014	estimated_generation_note_2015	estimated_generation_note_2016	estimated_generation_note_2017	geometry
0	AFG	Afghanistan	Kajaki Hydroelectric Power Plant Afghanistan	GEODB0040538	33.0	32.3220	65.1190	Hydro	NaN	NaN	...	162.90	97.39	137.76	119.50	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	POINT (65.11900 32.32200)
1	AFG	Afghanistan	Kandahar DOG	WKS0070144	10.0	31.6700	65.7950	Solar	NaN	NaN	...	17.48	18.25	17.70	18.29	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	POINT (65.79500 31.67000)
2	AFG	Afghanistan	Kandahar JOL	WKS0071196	10.0	31.6230	65.7920	Solar	NaN	NaN	...	17.58	19.10	17.62	18.72	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	POINT (65.79200 31.62300)
3	AFG	Afghanistan	Mahipar Hydroelectric Power Plant Afghanistan	GEODB0040541	66.0	34.5560	69.4787	Hydro	NaN	NaN	...	203.55	146.90	230.18	174.91	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	POINT (69.47870 34.55600)
4	AFG	Afghanistan	Naghlu Dam Hydroelectric Power Plant Afghanistan	GEODB0040534	100.0	34.6410	69.7170	Hydro	NaN	NaN	...	357.22	270.99	395.38	350.80	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	POINT (69.71700 34.64100)
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
34931	ZMB	Zambia	Ndola	WRI1022386	50.0	-12.9667	28.6333	Oil	NaN	NaN	...	NaN	NaN	NaN	183.79	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	CAPACITY-FACTOR-V1	POINT (28.63330 -12.96670)
34932	ZMB	Zambia	Nkana	WRI1022384	20.0	-12.8167	28.2000	Oil	NaN	NaN	...	NaN	NaN	NaN	73.51	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	CAPACITY-FACTOR-V1	POINT (28.20000 -12.81670)
34933	ZMB	Zambia	Victoria Falls	WRI1022380	108.0	-17.9167	25.8500	Hydro	NaN	NaN	...	575.78	548.94	579.90	578.32	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	POINT (25.85000 -17.91670)
34934	ZWE	Zimbabwe	Hwange Coal Power Plant Zimbabwe	GEODB0040404	920.0	-18.3835	26.4700	Coal	NaN	NaN	...	NaN	NaN	NaN	2785.10	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	CAPACITY-FACTOR-V1	POINT (26.47000 -18.38350)
34935	ZWE	Zimbabwe	Kariba Dam South Hydroelectric Power Station Z...	GEODB0003803	750.0	-16.5222	28.7619	Hydro	NaN	NaN	...	3999.75	3745.02	3428.65	3960.24	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	POINT (28.76190 -16.52220)

	country	country_long	name	gppd_idnr	capacity_mw	latitude	longitude	primary_fuel	other_fuel1	other_fuel2	...	estimated_generation_gwh_2014	estimated_generation_gwh_2015	estimated_generation_gwh_2016	estimated_generation_gwh_2017	estimated_generation_note_2013	estimated_generation_note_2014	estimated_generation_note_2015	estimated_generation_note_2016	estimated_generation_note_2017	geometry
17685	NLD	Netherlands	Ameland	WKS0066277	6.0	53.4503	5.6734	Solar	NaN	NaN	...	8.95	8.98	9.00	8.99	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	POINT (5.67340 53.45030)
17686	NLD	Netherlands	Andijk	WKS0072735	15.0	52.7370	5.2290	Solar	NaN	NaN	...	22.77	22.51	22.53	22.52	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	POINT (5.22900 52.73700)
17687	NLD	Netherlands	Beabuorren	WRI1005281	10.4	53.0417	5.5014	Wind	NaN	NaN	...	22.29	23.53	24.76	22.35	WIND-V1	WIND-V1	WIND-V1	WIND-V1	WIND-V1	POINT (5.50140 53.04170)
17688	NLD	Netherlands	Borssele Nuclear Power Station	WRI1019300	485.0	51.4312	3.7174	Nuclear	NaN	NaN	...	NaN	NaN	NaN	3203.82	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	CAPACITY-FACTOR-V1	POINT (3.71740 51.43120)
17689	NLD	Netherlands	Budel-Dorplein	WKS0071714	43.9	51.2360	5.6150	Solar	NaN	NaN	...	62.58	62.52	63.02	62.02	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	POINT (5.61500 51.23600)
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
17751	NLD	Netherlands	Westereems	WRI1005474	213.3	53.4375	6.8613	Wind	NaN	NaN	...	706.12	684.72	652.30	663.99	WIND-V1	WIND-V1	WIND-V1	WIND-V1	WIND-V1	POINT (6.86130 53.43750)
17752	NLD	Netherlands	Weststad Oosterhout	WRI1005476	15.0	51.6450	4.8597	Wind	NaN	NaN	...	36.62	36.68	33.10	35.03	WIND-V1	WIND-V1	WIND-V1	WIND-V1	WIND-V1	POINT (4.85970 51.64500)
17753	NLD	Netherlands	Willem Anna Polder	WRI1005477	9.0	51.6392	4.1292	Wind	NaN	NaN	...	17.80	18.96	17.17	18.51	WIND-V1	WIND-V1	WIND-V1	WIND-V1	WIND-V1	POINT (4.12920 51.63920)
17754	NLD	Netherlands	Zierikzee	WKS0073210	14.1	51.6400	3.9900	Solar	NaN	NaN	...	21.57	21.39	21.63	21.44	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	POINT (3.99000 51.64000)
17755	NLD	Netherlands	flushing	WRI1005431	2.3	51.4425	3.5736	Wind	NaN	NaN	...	NaN	NaN	5.88	7.70	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	WIND-V1	WIND-V1	POINT (3.57360 51.44250)

	NL	1	2	3	4	5	6
time
2013-01-01 00:00:00	12400.581906	1223.730363	2619.697268	2614.012629	3021.780489	2921.361158	0.0
2013-01-01 01:00:00	12139.240389	1197.940319	2564.487305	2558.922470	2958.096647	2859.793647	0.0
2013-01-01 02:00:00	12109.418967	1194.997443	2558.187351	2552.636187	2950.829747	2852.768239	0.0

Splitting EEZ and Netherlands into regions¶

II - Land Eligibility Analysis¶

Onwind¶

Offshore¶

Solar¶

III - Furthers parts¶

Onshore Capacity Factor¶

Offshore Capacity Factor¶

Solar Capacity Factor¶

PyPSA¶

Population¶

Costs¶

Pypsa Network¶

Buses creation¶

Carriers¶

Generators¶

Links and Storage¶

Optimization¶

Scenario No CO2 Emissions¶

Outputs of interest¶

Shadow Price¶

Curtailment¶

Sensitivity Analysis¶

Nuclear cost 5000€/kW¶

Nuclear cost 7500€/kW¶

Nuclear cost 10000€/kW¶

Sensitivity Analysis Results¶

	MRGID	GEONAME	MRGID_TER1	POL_TYPE	MRGID_SOV1	TERRITORY1	ISO_TER1	SOVEREIGN1	...	ISO_SOV1	ISO_SOV2	ISO_SOV3	UN_SOV1	UN_SOV2	UN_SOV3	UN_TER1	UN_TER2	UN_TER3	geometry
61	26519.0	Aruban Exclusive Economic Zone	8652.0	200NM	15.0	Aruba	ABW	Netherlands	...	NLD	None	None	528	NaN	NaN	533.0	NaN	NaN	MULTIPOLYGON (((-68.87017 15.24185, -68.88133 ...
66	26520.0	Bonaire Exclusive Economic Zone	15074.0	200NM	15.0	Bonaire	BES	Netherlands	...	NLD	None	None	528	NaN	NaN	535.0	NaN	NaN	MULTIPOLYGON (((-68.40522 11.66667, -68.42272 ...
74	26517.0	Curaçaoan Exclusive Economic Zone	19143.0	200NM	15.0	Curaçao	CUW	Netherlands	...	NLD	None	None	528	NaN	NaN	531.0	NaN	NaN	MULTIPOLYGON (((-68.86222 15.24111, -68.55902 ...
92	26518.0	Saba Exclusive Economic Zone	19095.0	200NM	15.0	Saba	BES	Netherlands	...	NLD	None	None	528	NaN	NaN	535.0	NaN	NaN	MULTIPOLYGON (((-63.28635 17.87370, -63.26729 ...
94	26526.0	Sint-Eustatius Exclusive Economic Zone	26525.0	200NM	15.0	Sint-Eustatius	BES	Netherlands	...	NLD	None	None	528	NaN	NaN	535.0	NaN	NaN	MULTIPOLYGON (((-62.76955 17.64394, -62.80457 ...
100	21803.0	Sint-Maarten Exclusive Economic Zone	21802.0	200NM	15.0	Sint-Maarten	SXM	Netherlands	...	NLD	None	None	528	NaN	NaN	534.0	NaN	NaN	MULTIPOLYGON (((-63.01457 18.05186, -63.01409 ...
216	5668.0	Dutch Exclusive Economic Zone	15.0	200NM	15.0	Netherlands	NLD	Netherlands	...	NLD	None	None	528	NaN	NaN	528.0	NaN	NaN	MULTIPOLYGON (((7.21028 53.24417, 7.20750 53.2...

	GID_0	GID_1	COUNTRY	NAME_1	geometry
0	NLD	NLD.1_1	Netherlands	Drenthe	MULTIPOLYGON (((6.52366 52.64551, 6.52366 52.6...
1	NLD	NLD.2_1	Netherlands	Flevoland	MULTIPOLYGON (((5.30388 52.31247, 5.28622 52.3...
2	NLD	NLD.3_1	Netherlands	Fryslân	MULTIPOLYGON (((6.23919 52.91320, 6.23615 52.9...
3	NLD	NLD.4_1	Netherlands	Gelderland	MULTIPOLYGON (((5.13576 51.73854, 5.13544 51.7...
4	NLD	NLD.5_1	Netherlands	Groningen	MULTIPOLYGON (((6.91710 53.01160, 6.91460 53.0...
5	NLD	NLD.6_1	Netherlands	IJsselmeer	MULTIPOLYGON (((5.26216 52.43537, 5.23710 52.4...
6	NLD	NLD.7_1	Netherlands	Limburg	MULTIPOLYGON (((6.01553 50.74829, 6.01326 50.7...
7	NLD	NLD.8_1	Netherlands	Noord-Brabant	MULTIPOLYGON (((5.50305 51.29579, 5.49979 51.2...
8	NLD	NLD.9_1	Netherlands	Noord-Holland	MULTIPOLYGON (((4.83937 52.23358, 4.83403 52.2...
9	NLD	NLD.10_1	Netherlands	Overijssel	MULTIPOLYGON (((6.87309 52.13272, 6.87284 52.1...
10	NLD	NLD.11_1	Netherlands	Utrecht	MULTIPOLYGON (((5.15190 51.93248, 5.15337 51.9...
11	NLD	NLD.12_1	Netherlands	Zeeland	MULTIPOLYGON (((3.67419 51.29042, 3.66687 51.2...
12	NLD	NLD.13_1	Netherlands	Zeeuwse meren	MULTIPOLYGON (((3.86646 51.53969, 3.85926 51.5...
13	NLD	NLD.14_1	Netherlands	NA	MULTIPOLYGON (((4.36117 51.66149, 4.33181 51.6...

	geometry	GID_0	GID_1	COUNTRY	NAME_1
region
1	MULTIPOLYGON (((6.10234 52.84489, 6.10211 52.8...	NLD	NLD.1_1	Netherlands	Drenthe
2	MULTIPOLYGON (((5.13285 51.73853, 5.13139 51.7...	NLD	NLD.2_1	Netherlands	Flevoland
3	POLYGON ((6.01326 50.75231, 6.01215 50.75426, ...	NLD	NLD.7_1	Netherlands	Limburg
4	MULTIPOLYGON (((5.14756 51.92529, 5.14031 51.9...	NLD	NLD.9_1	Netherlands	Noord-Holland
5	MULTIPOLYGON (((3.66083 51.29253, 3.65546 51.2...	NLD	NLD.12_1	Netherlands	Zeeland
6	MULTIPOLYGON (((3.85372 51.53911, 3.84910 51.5...	NLD	NLD.6_1	Netherlands	IJsselmeer

	geometry	GID_0	GID_1	COUNTRY	NAME_1	MRGID	GEONAME	MRGID_TER1	POL_TYPE	MRGID_SOV1	...	AREA_KM2	ISO_SOV1	ISO_SOV2	ISO_SOV3	UN_SOV1	UN_SOV2	UN_SOV3	UN_TER1	UN_TER2	UN_TER3
region
1	MULTIPOLYGON (((5.27139 53.25236, 5.27139 53.2...	NLD	NLD.1_1	Netherlands	Drenthe	NaN	None	NaN	None	NaN	...	NaN	None	None	None	NaN	NaN	NaN	NaN	NaN	NaN
2	MULTIPOLYGON (((5.12990 51.74391, 5.13257 51.7...	NLD	NLD.2_1	Netherlands	Flevoland	NaN	None	NaN	None	NaN	...	NaN	None	None	None	NaN	NaN	NaN	NaN	NaN	NaN
3	POLYGON ((6.01326 50.75231, 6.01215 50.75426, ...	NLD	NLD.7_1	Netherlands	Limburg	NaN	None	NaN	None	NaN	...	NaN	None	None	None	NaN	NaN	NaN	NaN	NaN	NaN
4	MULTIPOLYGON (((4.99422 52.36170, 4.99569 52.3...	NLD	NLD.9_1	Netherlands	Noord-Holland	NaN	None	NaN	None	NaN	...	NaN	None	None	None	NaN	NaN	NaN	NaN	NaN	NaN
5	MULTIPOLYGON (((3.64316 51.29162, 3.63084 51.2...	NLD	NLD.12_1	Netherlands	Zeeland	NaN	None	NaN	None	NaN	...	NaN	None	None	None	NaN	NaN	NaN	NaN	NaN	NaN
6	MULTIPOLYGON (((4.22318 51.51344, 4.22317 51.5...	NLD	NLD.6_1	Netherlands	IJsselmeer	5668.0	Dutch Exclusive Economic Zone	15.0	200NM	15.0	...	64328.0	NLD	None	None	528.0	NaN	NaN	528.0	NaN	NaN

region	6
time
2013-01-01 00:00:00	0.728403
2013-01-01 01:00:00	0.657270
2013-01-01 02:00:00	0.641306
2013-01-01 03:00:00	0.574383
2013-01-01 04:00:00	0.443536
...	...
2013-12-31 19:00:00	0.999919
2013-12-31 20:00:00	0.992884
2013-12-31 21:00:00	0.974395
2013-12-31 22:00:00	0.960512
2013-12-31 23:00:00	0.922488

primary_fuel	Coal	Gas	Nuclear	Solar	Wind
region
1.0	1600.0	4306.0	0.0	62.8	320.18
2.0	590.0	349.0	0.0	0.0	49.10
3.0	0.0	0.0	0.0	67.9	35.80
4.0	685.0	1025.0	0.0	30.0	38.73
5.0	2140.0	1753.0	485.0	80.3	174.28
6.0	0.0	1879.0	0.0	0.0	1178.10

	population	region	ratio
Drenthe	502051	1	0.028187
Flevoland	444701	2	0.024967
Fryslân	659551	1	0.037030
Gelderland	2133708	2	0.119795
Groningen	596075	1	0.033466
IJsselmeer	0	6	0.000000
Limburg	1128367	3	0.063351
Noord-Brabant	2626210	3	0.147446
Noord-Holland	2952622	4	0.165772
Overijssel	1184333	2	0.066493
Utrecht	1387643	4	0.077908
Zeeland	391124	5	0.021959
Zeeuwse meren	0	6	0.000000
NA	3804906	5	0.213623
Dutch Exclusive Economic Zone	0	6	0.000000

	ratio	population
region
1	0.098683	1757677
2	0.211256	3762742
3	0.210798	3754577
4	0.243681	4340265
5	0.235583	4196030
6	0.000000	0

		value	unit	source	further description	currency_year
technology	parameter
Ammonia cracker	FOM	4.300000e+00	%/year	Ishimoto et al. (2020): 10.1016/j.ijhydene.202...	Estimated based on Labour cost rate, Maintenan...	2015.0
	ammonia-input	1.460000e+00	MWh_NH3/MWh_H2	ENGIE et al (2020): Ammonia to Green Hydrogen ...	Assuming a integrated 200t/d cracking and puri...	NaN
	investment	1.123945e+06	EUR/MW_H2	Ishimoto et al. (2020): 10.1016/j.ijhydene.202...	Calculated. For a small (200 t_NH3/d input) fa...	2015.0
	lifetime	2.500000e+01	years	Ishimoto et al. (2020): 10.1016/j.ijhydene.202...	NaN	2015.0
BEV Bus city	FOM	2.000000e-04	%/year	Danish Energy Agency, inputs/data_sheets_for_c...	BEV B1	2022.0
...	...	...	...	...	...	...
waste CHP CC	efficiency-heat	7.627000e-01	per unit	Danish Energy Agency, technology_data_for_el_a...	08 WtE CHP, Large, 50 degree: Heat efficiency...	2015.0
	investment	8.829851e+03	EUR/kW_e	Danish Energy Agency, technology_data_for_el_a...	08 WtE CHP, Large, 50 degree: Nominal investm...	2015.0
	lifetime	2.500000e+01	years	Danish Energy Agency, technology_data_for_el_a...	08 WtE CHP, Large, 50 degree: Technical lifetime	2015.0
water tank charger	efficiency	8.367000e-01	per unit	Danish Energy Agency, technology_data_catalogu...	: efficiency from sqr(Round trip efficiency)	2015.0
water tank discharger	efficiency	8.367000e-01	per unit	Danish Energy Agency, technology_data_catalogu...	: efficiency from sqr(Round trip efficiency)	2015.0

parameter	C in fuel	C stored	CO2 intensity	CO2 stored	FOM	Motor size	VOM	ammonia-input	c_b	c_v	...	lifetime	lohc-input	methane-input	methanol-input	min_fill_level	naphtha-input	nitrogen-input	ore-input	p_nom_ratio	pelletizing cost
technology
Ammonia cracker	NaN	NaN	0.0	NaN	4.3000	NaN	0.0000	1.46	NaN	NaN	...	25.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
BEV Bus city	NaN	NaN	0.0	NaN	0.0002	320.6897	0.0952	NaN	NaN	NaN	...	12.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
BEV Coach	NaN	NaN	0.0	NaN	0.0001	298.2759	0.0952	NaN	NaN	NaN	...	12.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
BEV Truck Semi-Trailer max 50 tons	NaN	NaN	0.0	NaN	0.0002	468.9655	0.0952	NaN	NaN	NaN	...	10.5	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
BEV Truck Solo max 26 tons	NaN	NaN	0.0	NaN	0.0002	353.4483	0.0952	NaN	NaN	NaN	...	13.8	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
uranium	NaN	NaN	0.0	NaN	0.0000	NaN	0.0000	NaN	NaN	NaN	...	25.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
waste CHP	NaN	NaN	0.0	NaN	2.3789	NaN	28.4644	NaN	0.2872	1.0	...	25.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
waste CHP CC	NaN	NaN	0.0	NaN	2.3789	NaN	28.4644	NaN	0.2872	1.0	...	25.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
water tank charger	NaN	NaN	0.0	NaN	0.0000	NaN	0.0000	NaN	NaN	NaN	...	25.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
water tank discharger	NaN	NaN	0.0	NaN	0.0000	NaN	0.0000	NaN	NaN	NaN	...	25.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

	geometry	x	y
region
1	POINT (6.30770 53.03945)	6.307700	53.039446
2	POINT (6.29240 52.29567)	6.292401	52.295671
3	POINT (5.83110 51.27609)	5.831101	51.276089
4	POINT (4.80150 52.41195)	4.801503	52.411945
5	POINT (4.41346 51.97036)	4.413465	51.970355
6	POINT (4.13744 53.40210)	4.137437	53.402103

attribute	bus	control	p_nom	p_nom_mod	p_nom_extendable	p_nom_min	p_nom_max	p_min_pu	p_max_pu	...	min_up_time	min_down_time	up_time_before	down_time_before	ramp_limit_up	ramp_limit_down	ramp_limit_start_up	ramp_limit_shut_down	weight	p_nom_opt
Generator
coal 1	1	PQ	1600.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
gas 1	1	PQ	4306.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
nuclear 1	1	PQ	0.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
coal 2	2	PQ	590.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
gas 2	2	PQ	349.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
nuclear 2	2	PQ	0.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
coal 3	3	PQ	0.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
gas 3	3	PQ	0.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
nuclear 3	3	PQ	0.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
coal 4	4	PQ	685.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
gas 4	4	PQ	1025.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
nuclear 4	4	PQ	0.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
coal 5	5	PQ	2140.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
gas 5	5	PQ	1753.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
nuclear 5	5	PQ	485.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
coal 6	6	PQ	0.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
gas 6	6	PQ	1879.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
nuclear 6	6	PQ	0.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0

	geometry	x	y
region
1	POINT (4073443.665 3331958.052)	6.307700	53.039446
2	POINT (4068174.551 3249337.655)	6.292401	52.295671
3	POINT (4030273.524 3137753.671)	5.831101	51.276089
4	POINT (3967606.828 3268464.915)	4.801503	52.411945
5	POINT (3937508.981 3221406.627)	4.413465	51.970355
6	POINT (3931522.763 3381771.783)	4.137437	53.402103

Bus	1	2	3	4	5	6
snapshot
2013-01-01 00:00:00	24.568	1.508	941.472490	550.600443	32.213321	0.0
2013-01-01 01:00:00	24.568	1.508	941.472490	550.600443	32.213321	0.0
2013-01-01 02:00:00	24.568	1.508	941.472490	603.712006	32.213321	0.0
2013-01-01 03:00:00	24.568	1.508	941.472490	603.712006	32.213321	0.0
2013-01-01 04:00:00	24.568	1.508	941.472490	603.712006	32.213321	0.0
...	...	...	...	...	...	...
2013-12-31 19:00:00	24.568	1.508	1032.288029	603.712006	32.213321	0.0
2013-12-31 20:00:00	24.568	1.508	1032.288029	603.712006	32.213321	0.0
2013-12-31 21:00:00	24.568	1.508	1032.288029	603.712006	32.213321	0.0
2013-12-31 22:00:00	24.568	1.508	1032.288029	603.712006	32.213321	0.0
2013-12-31 23:00:00	24.568	1.508	1032.288029	603.712006	32.213321	0.0

	value	unit	source	further description	currency_year
parameter
CO2 intensity	0.198	tCO2/MWh_th	Stoichiometric calculation with 50 GJ/t CH4	NaN	NaN
fuel	24.568	EUR/MWh_th	DIW (2013): Current and propsective costs of e...	Based on IEA 2011 data.	2010.0

parameter	C in fuel	C stored	CO2 intensity	CO2 stored	FOM	Motor size	VOM	ammonia-input	c_b	c_v	...	methane-input	methanol-input	min_fill_level	naphtha-input	nitrogen-input	ore-input	p_nom_ratio	pelletizing cost	marginal_cost	capital_cost
technology
nuclear	NaN	NaN	0.0	NaN	1.2700	NaN	3.5464	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	14.013271	753784.215297
offwind	NaN	NaN	0.0	NaN	2.3741	NaN	0.0212	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.021200	184567.393550
onwind	NaN	NaN	0.0	NaN	1.2347	NaN	1.5080	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	1.508000	107441.457021
solar	NaN	NaN	0.0	NaN	1.7275	NaN	0.0106	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.010600	63114.923985